library(dplyr)
library(caret)
library(ggplot2)
library(lattice)
library(class)
library(corrplot)
library(pROC)
library(caret)
library(lattice)
library(foreach)
library(iterators)
library(parallel)
library(doMC)
library(e1071)
library(ranger)
library(rpart)
library(rpart.plot)
library(NbClust)
library(psych)
library(factoextra)
library(ade4)
library(cluster)
library(lift)
library(dplyr)
library(readr)
library(mice)
library(stringr)
library(ggplot2)
library(tidyverse)
library(gamlss)
library(MASS)
library(doMC)
library(iterators)
library(gam)
library(stats)
library(splines)
library(foreach)
library(gamlss.data)
library(gamlss.dist)
library(NB)
library(gridExtra)
library(gmodels)
library(Hmisc)
library(ggthemes)
library(ellipse)
library(ggcorrplot)
library(corrplot)
library(leaps)
library(PerformanceAnalytics)
library(GGally)
library(psych)
library(carData)
library(car)
library(lmtest)
library(olsrr)
library(performance)
library(mice)
library(see)
library(lme4)
library(patchwork)

Introducción

En esta sección se comenzarán a aplicar los algoritmos de aprendizaje automático (o machine learning) vistos en clase. Utilizaremos los siguientes algoritmos que podemos clasificar en:

  • Algoritmos no supervisados: Parten de un conjunto de datos sin etiquetar y aplican inducción partiendo de ejemplos para predecir en nuevos datos.

    • Análisis de componentes principales (reducción de la dimensionalidad)
    • k-medias. Análisis cluster no jerárquico o de conglomerados (clustering)
    • Análisis cluster jerárquico o de conglomerados (clustering)
  • Algoritmos supervisados: Utilizan un conjunto de datos etiquetados y con ellas intentan encontrar relaciones entre las características y las etiquetas para predecir en nuevos datos.

    • Regresión logística
    • Los k-vecinos más cercanos (k-NN: The k-nearest neighbours) (Basado en la información)
    • Árboles de decisión (Basado en la información)
    • Bosques aleatorios (Basado en la información)
    • Máquinas de vectores de soporte (SVM: Support Vector Machines) (Basado en el error)

Se carga, limpia, visualiza y secciona el conjunto de datos:

data <- read.csv("train-ML.csv", na = c("","NA","NULL",NULL,"  ","/n" ))
head(data)
dataTest <- read.csv("test-ML.csv", na = c("","NA","NULL",NULL,"  ","/n" ))
head(dataTest)
library(dplyr)
data %>% dplyr::select(-Price.mod2) -> data
data %>% dplyr::select(-X) -> data
head(data)
data %>%
  mutate(Fuel_Type = as.factor(Fuel_Type)) %>%
  mutate(Transmission = as.factor(Transmission)) %>%
  mutate(Location = as.factor(Location)) -> data
data %>% mutate(Owner_Type = factor(Owner_Type, levels=c("First", "Second", "Third", "Fourth & Above"))) -> data

dataTest %>%
  mutate(Fuel_Type = as.factor(Fuel_Type)) %>%
  mutate(Transmission = as.factor(Transmission)) %>%
  mutate(Location = as.factor(Location)) -> dataTest

dataTest %>%
  mutate(Name = as.factor(Name)) %>%
  mutate(Make = as.factor(Make)) %>%
  mutate(Seats = as.factor(Seats))%>%
  mutate(Gama = as.factor(Gama))-> dataTest

data %>%
  mutate(Name = as.factor(Name)) %>%
  mutate(Make = as.factor(Make)) %>%
  mutate(Seats = as.factor(Seats))-> data

dataTest %>% mutate(Owner_Type = factor(Owner_Type, levels=c("First", "Second", "Third", "Fourth & Above"))) -> dataTest

Análisis no supervisado

En primer lugar, comenzaremos viendo los algoritmos de aprendizaje no supervisado.

Análisis de componentes principales

El análisis de componentes principales (en inglés, PCA) es una técnica utilizada para describir un conjunto de datos en términos de nuevas variables denominadas componentes no correlacionadas. Estas nuevas componentes se construyen a partir de las variables existentes, eso sí, debemos asegurarnos de que las variables utilizadas en PCA sean variables cuantitativas (no podemos usar variables cualitativas ni categóricas). Con esta técnica se pretende reducir la dimensionalidad del problema en cuestión.

str(data)
## 'data.frame':    4711 obs. of  15 variables:
##  $ Name             : Factor w/ 1636 levels "Ambassador Classic Nova Diesel",..: 37 72 501 1103 443 221 1073 1248 290 121 ...
##  $ Location         : Factor w/ 11 levels "Ahmedabad","Bangalore",..: 10 10 3 8 11 5 10 4 10 3 ...
##  $ Year             : int  2015 2014 2016 2009 2005 2016 2016 2017 2011 2012 ...
##  $ Kilometers_Driven: int  48000 18600 18000 80464 123200 46727 17000 30090 32000 115000 ...
##  $ Fuel_Type        : Factor w/ 4 levels "CNG","Diesel",..: 2 2 4 2 4 2 4 2 4 2 ...
##  $ Transmission     : Factor w/ 2 levels "Automatic","Manual": 1 1 1 1 2 2 1 1 2 1 ...
##  $ Owner_Type       : Factor w/ 4 levels "First","Second",..: 1 2 1 2 2 1 1 1 2 3 ...
##  $ Engine           : int  1968 1995 1197 2987 1495 1498 1595 1461 1196 1995 ...
##  $ Power            : num  174 190 82 198 94 ...
##  $ Seats            : Factor w/ 8 levels "2","4","5","6",..: 3 3 3 3 3 3 3 3 3 3 ...
##  $ Price            : num  18.25 21 5.4 9.29 1 ...
##  $ Make             : Factor w/ 28 levels "Ambassador","Audi",..: 2 4 11 18 11 9 18 23 9 4 ...
##  $ Gama             : chr  "Gama alta" "Gama alta" "Gama media" "Gama alta" ...
##  $ kmpl             : num  15.7 22.7 18.9 11 13.2 ...
##  $ kmpkg            : num  0 0 0 0 0 0 0 0 0 0 ...
data
corrplot(cor(data[, c(3,4, 8, 9, 14, 15)]), method = "ellipse") 

corPlot(data[, c(3,4, 8, 9, 14, 15)], cex = 1.2, main = "Matriz de correlación")

corrplot(cor(data[, c(3,4, 8, 9, 14, 15)]),method = "circle",       order = "hclust",         hclust.method = "ward.D",
         addrect =2,rect.col = 3,rect.lwd = 3)  

cortest(cor(data[, c(3,4, 8, 9, 14, 15)]))
## Tests of correlation matrices 
## Call:cortest(R1 = cor(data[, c(3, 4, 8, 9, 14, 15)]))
##  Chi Square value 250.68  with df =  15   with probability < 9e-45

Tenemos evidencias para decir que las correlaciones son distintas de 0.

pca1 <- prcomp(data[, c(3,4, 8, 9, 14, 15)])
plot(pca1)

summary(pca1)
## Importance of components:
##                           PC1       PC2   PC3   PC4   PC5   PC6
## Standard deviation     101058 600.56073 27.03 4.409 3.046 1.968
## Proportion of Variance      1   0.00004  0.00 0.000 0.000 0.000
## Cumulative Proportion       1   1.00000  1.00 1.000 1.000 1.000

En nuestro conjunto de datos inicial, todas nuestras variables eran cuantitativas (menos la variable respuesta, que no utilizamos en aprendizaje no supervisado), sin embargo, las variables categorizadas que hemos creado no lo son. Así que haremos el análisis de componentes principales escalando las variables: Kilometers_Driven, Power, kmpl y kmpg.

pca2 <- prcomp(data[, c(3,4, 8, 9, 14, 15)], scale=T)

pca2
## Standard deviations (1, .., p=6):
## [1] 1.5099667 1.1444561 1.0599891 0.9273837 0.5500460 0.3522111
## 
## Rotation (n x k) = (6 x 6):
##                           PC1        PC2          PC3        PC4         PC5
## Year              -0.13713864 -0.3782024  0.596223411 0.63439205  0.28287926
## Kilometers_Driven  0.10115280  0.1465697 -0.668955170 0.71929106  0.04377908
## Engine             0.61114224 -0.2395341  0.007784461 0.01713324 -0.11602069
## Power              0.58101040 -0.3074207  0.085329519 0.04975366 -0.38891068
## kmpl              -0.50571711 -0.4201464 -0.135355634 0.08348424 -0.71719694
## kmpkg              0.06436916  0.7120907  0.413948973 0.26538370 -0.48885474
##                           PC6
## Year              -0.01412815
## Kilometers_Driven  0.03857703
## Engine            -0.74519368
## Power              0.63789605
## kmpl              -0.16752599
## kmpkg             -0.08956703

Aquí arriba, podemos ver los diferentes pesos que otorga el análisis de componentes principales a cada una de las variables iniciales escaladas. Por ejemplo: en la primera componente principal (PC1), vemos que sobre todo se enfrentan las variables Engine y Power contra kmpl; en la segunda componente principal (PC2), vemos que se enfrenta la variable kmpg contra kmpl, Power,Year y Engine; en la tercera componente principal (PC3) se enfrentan los kilómetros que lleva recorridos el coche contra el año de fabricación y la variable kmpkg.

summary(pca2)
## Importance of components:
##                         PC1    PC2    PC3    PC4     PC5     PC6
## Standard deviation     1.51 1.1445 1.0600 0.9274 0.55005 0.35221
## Proportion of Variance 0.38 0.2183 0.1873 0.1433 0.05043 0.02068
## Cumulative Proportion  0.38 0.5983 0.7856 0.9289 0.97932 1.00000

La inercia de las primeras dimensiones muestra si existen relaciones fuertes entre las variables y sugiere el número de dimensiones que se deben estudiar.

Las dos primeras dimensiones de análisis expresan el 59,83% de la inercia total del conjunto de datos; eso quiere decir que el 59.83% de los individuos (o variables) nublan la variabilidad total que es explicada por el plano. Este porcentaje es relativamente alto y, por lo tanto, el primer plano representa bien la variabilidad de los datos. Este valor es muy superior al valor de referencia que equivale al 34,95%, por lo que la variabilidad explicada por este plano es muy significativa (el valor de referencia es el cuantil 0,95 de la distribución de porcentajes de inercia obtenida simulando 1689 tablas de datos de tamaño equivalente sobre la base de una distribución normal).

A partir de estas observaciones, conviene interpretar también las dimensiones mayores o iguales a la tercera.

Sin embargo, aquí resulta difícil ver algo claro e intuitivo, así que haremos un pequeño resumen y un gráfico multivariante para mostrar la información más relevante del PCA. Se estudia el plano 1:2.

PC1= pca2[[2]][,1]
PC2= pca2[[2]][,2]
PC3= pca2[[2]][,3]
PC4= pca2[[2]][,4]

componentes_princ <- cbind(PC1,PC2,PC3,PC4)
componentes_princ
##                           PC1        PC2          PC3        PC4
## Year              -0.13713864 -0.3782024  0.596223411 0.63439205
## Kilometers_Driven  0.10115280  0.1465697 -0.668955170 0.71929106
## Engine             0.61114224 -0.2395341  0.007784461 0.01713324
## Power              0.58101040 -0.3074207  0.085329519 0.04975366
## kmpl              -0.50571711 -0.4201464 -0.135355634 0.08348424
## kmpkg              0.06436916  0.7120907  0.413948973 0.26538370
plot(pca2)

biplot(pca2)

s.corcircle(componentes_princ[,-3], sub="PC1 Y PC2")

En el gráfico enfrentamos la primera y la segunda componente principal, y vemos como influyen cada una de las variables en los coches. Por ejemplo, el coche 1712, debe tener unos valores muy altos de Power y Engine que son las variables que “más tiran hacia la derecha”, y el coche 745, según el biplot debe tener un valor muy alto de kilometers_driven, que es la variable que “más tira en esa dirección”. Si recordamos lo analizado previamente, en la sección del análisis exploratorio de datos, este coche ya destacó por tener valores un tanto diferentes a los del resto por su desmesurado valor de kilómetros recorridos. Comprobamos que la información que nos proporciona el gráfico es totalmente coherente con lo que obtuvimos en el EDA.

La dimensión 1 opone individuos caracterizados por una coordenada fuertemente positiva en el eje (a la derecha del gráfico) a individuos caracterizados por una coordenada fuertemente negativa en el eje (a la izquierda del gráfico).

El grupo 1 (caracterizado por una coordenada positiva en el eje) comparte:

  • valores altos para la variable kmpkg.
  • valores bajos para las variables kmpl, Potencia y Motor (las variables se ordenan de la más débil).

El grupo 2 (caracterizado por una coordenada positiva en el eje) comparte:

  • valores altos para las variables Motor, Potencia y Kilómetros_Recorridos (las variables se ordenan de la más fuerte).
  • valores bajos para las variables kmpl, Year y kmpkg (las variables se ordenan de la más débil).

El grupo 3 (caracterizado por una coordenada negativa en el eje) comparte:

  • valores altos para las variables kmpl y Año (las variables se ordenan de la más fuerte).
  • valores bajos para las variables Motor, Potencia, Kilómetros_Conducidos y kmpkg (las variables se ordenan de menor a mayor).

La dimensión 2 opone individuos caracterizados por una coordenada fuertemente positiva en el eje (en la parte superior del gráfico) a individuos caracterizados por una coordenada fuertemente negativa en el eje (en la parte inferior del gráfico).

El grupo 1 (caracterizado por una coordenada positiva en el eje) comparte:

  • valores altos para las variables kmpl y Año (las variables se ordenan de la más fuerte).
  • valores bajos para las variables Motor, Potencia, Kilómetros_Conducidos y kmpkg (las variables se ordenan de menor a mayor).

El grupo 2 (caracterizado por una coordenada positiva en el eje) comparte:

  • valores altos para las variables Motor, Potencia y Kilómetros_Recorridos (las variables se ordenan de la más fuerte).
  • valores bajos para las variables kmpl, Year y kmpkg (las variables se ordenan de la más débil).

El grupo 3 (caracterizado por una coordenada negativa en el eje) comparte:

  • valores altos para la variable kmpkg.
  • valores bajos para las variables kmpl, Potencia y Motor (las variables se ordenan de la más débil).
fviz_pca_var(pca2,axes = c(1,2), col.var = "cos2", alpha.var = "contrib" ) + theme_grey()

fviz_pca_var(pca2,axes = c(1,3), col.var = "cos2", alpha.var = "contrib" ) + theme_grey()

La suma de cos2 de una variable determinada sobre cada factor es 1. Esto significa que cada vector debería estar tocando el perímetro de la circunferencia unidad, pero no lo está haciendo ninguna prácticamente, ¿por qué?. Si observamos por ejemplo la variable Engine(al igual que Power), vemos que está muy cerca de tocar dicho perímetro, su proyección sobre las dimensiones 1 y 2 (componentes) indica su contribución a éstas, pero aún le falta algo de contribución que debe estar repartida por otra u otras dimensiones. Si está variable solo tuviese peso sobre las dos primeras dimensiones estaría tocando la circunferencia.

Podemos colorear las observaciones según alguna variable. Además podemos hacer que las variables que más contribuyen en este plano factorial, se resalten más que las que menos influencia tienen. También tenemos la posibilidad de dibujar elipses alrededor de cada grupo con un cierto nivel de confianza.

Como se aprecia en los gráficos anteriores, no tiene mucho sentido representar la segunda componente principal ya que no realiza un correcto enfrentamiento de variables y no nos aporta nada.

summary(pca2)
## Importance of components:
##                         PC1    PC2    PC3    PC4     PC5     PC6
## Standard deviation     1.51 1.1445 1.0600 0.9274 0.55005 0.35221
## Proportion of Variance 0.38 0.2183 0.1873 0.1433 0.05043 0.02068
## Cumulative Proportion  0.38 0.5983 0.7856 0.9289 0.97932 1.00000

Por otro lado, podemos decir que lo que más nos interesa de este resumen es la proporción de la varianza total que consigue explicar cada componente principal. Según el resumen que acabamos de mostrar arriba, vemos que la varianza total explicada no aumenta mucho a partir de la tercera o cuarta componente principal (y que con todas las componentes principales, evidentemente, la varianza explicada es el 100%). Para visualizar esto haremos un gráfico de barras:

screeplot(pca2, xlab="PCs")

Una estimación del número correcto de ejes a interpretar sugiere restringir el análisis a la descripción de los 3 primeros ejes. Estos ejes presentan una inercia superior a las obtenidas por el cuantil 0,95 de las distribuciones aleatorias (78,56% frente a 51,79%). Esta observación sugiere que solo estos ejes llevan una información real. En consecuencia, la descripción se situará en estos ejes.

scree(data[, c(3,4, 8, 9, 14, 15)],main ="Grafico_de_Sedimentacion")

El grafico de sedimentación nos muestra la cantidad óptima de componentes a tomar en el análisis, siendo los valores por encima de la linea de 1.0 los más aceptables.

fa.parallel(data[, c(3,4, 8, 9, 14, 15)],fa="pc")

## Parallel analysis suggests that the number of factors =  NA  and the number of components =  3

Según los resultados del análisis paralelo, el número de componentes deberá ser 3.

Se comprueba que con PCA no se consigue lo que se busca, ni PC2 ni PC3 nos sirven para realizar una correcta redimensión de los datos. Nos damos cuenta de que nuestro problema es bastante difícil de resolver, dado que es complicado ver algún tipo de separación o tendencia de los coches en función de las variables explicativas o incluso en cuanto al precio. Aún así, si nos fijamos, en los gráficos que enfrentan PC1 con PC2, PC1 con PC3 y PC2 con PC3, parece que los coches de gama media y gama alta están más dispersos que los de gama baja.

k-medias. Análisis cluster no jerárquico o de conglomerados (clustering)

El análisis cluster busca agrupar individuos u observaciones tratando de lograr la máxima homogeneidad en cada grupo o cluster y la mayor diferencia entre los grupos. Es decir, el clustering asigna individuos similares al mismo grupo, y observaciones muy distintas estarán en grupos diferentes. Nosotros usaremos el algoritmo de las k-medias que tiene como objetivo encontrar y agrupar en clases a las observaciones que tienen una alta similitud entre ellos. Esta similitud se define como la menor distancia entre características de cada observación. Cuanto más cerca estén los puntos de datos, más similares serán y más probabilidad habrá de que pertenezcan al mismo cluster. Para ello, primero debemos escoger una distancia y dado que nuestra variable respuesta está bastante bien balanceada, usaremos la distancia euclídea. Esta es la distancia que muchos métodos de R utilizan por defecto, pero debemos asegurarnos de que los datos que introducimos en el algoritmo están escalados (para que no tengan mayor importancia las variables con números más grandes en valor absoluto, por el mero de hecho de que puedan estar medidas en diferentes unidades, por ejemplo). Así y para trabajar sobre una variable respuesta binaria, transformamos la variable Gama.

data <- data[,-13]
data %>% 
  mutate(
    Gama = case_when( 
      data$Make=="Datsun" |data$Make=="Smart" |data$Make=="Tata" |data$Make=="Fiat" |data$Make=="Chevrolet" |data$Make=="Ambassador" |
      data$Make=="Skoda"|data$Make=="Renault" |data$Make=="Ford" |data$Make=="Honda"|data$Make=="Volkswagen" |data$Make=="Hyundai" |data$Make=="Nissan" |data$Make=="Maruti" ~ "Baja-media",
      data$Make=="Bentley"|data$Make=="Porsche" |data$Make=="Land Rover" |data$Make=="Jaguar" |data$Make=="Mini" |data$Make=="Mercedes-Benz" |data$Make=="Audi" |data$Make=="Bmw" |data$Make=="Jeep" |data$Make=="Volvo" |data$Make=="Isuzu" |data$Make=="Mitsubishi" |data$Make=="Toyota" |data$Make=="Force" |data$Make=="Mahindra"~ "Alta"
    )
  ) -> data

data$Gama <- as.factor(data$Gama)

#dataTest %>% dplyr::select(-Mileage) -> dataTest
dataTest %>% dplyr::select(-X) -> dataTest
head(dataTest)
dataTest%>%drop_na() -> dataTest
dataTest %>% 
    mutate(kmpl = ifelse(Fuel_Type=="Diesel" | Fuel_Type=="Petrol", Mileage, 0)) %>%
    mutate(kmpkg = ifelse(Fuel_Type=="CNG" | Fuel_Type=="LPG", Mileage, 0)) %>%
    dplyr::select(-Mileage) -> dataTest

dataTest %>% 
  mutate(
    Gama = case_when( 
      dataTest$Make=="Datsun" |dataTest$Make=="Smart" |dataTest$Make=="Tata" |dataTest$Make=="Fiat" |dataTest$Make=="Chevrolet" |dataTest$Make=="Ambassador" |
      dataTest$Make=="Skoda"|dataTest$Make=="Renault" |dataTest$Make=="Ford" |dataTest$Make=="Honda"|dataTest$Make=="Volkswagen" |dataTest$Make=="Hyundai" |dataTest$Make=="Nissan" |dataTest$Make=="Maruti" ~ "Baja-media",
      dataTest$Make=="Bentley"|dataTest$Make=="Porsche" |dataTest$Make=="Land Rover" |dataTest$Make=="Jaguar" |dataTest$Make=="Mini" |dataTest$Make=="Mercedes-Benz" |dataTest$Make=="Audi" |dataTest$Make=="Bmw" |dataTest$Make=="Jeep" |dataTest$Make=="Volvo" |dataTest$Make=="Isuzu" |dataTest$Make=="Mitsubishi" |dataTest$Make=="Toyota" |dataTest$Make=="Force" |dataTest$Make=="Mahindra"~ "Alta"
    )
  ) -> dataTest
TrainEscalado <- data %>% dplyr::select(Year, Kilometers_Driven, Engine, Power, Price, kmpl) %>% scale() %>% as.data.frame()

Como mínimo haremos dos grupos, es decir, buscaremos hacer 2 o más grupos, porque hacer un único cluster no tiene ningún sentido, ya que buscamos separar los coches en una característica que toma dos valores: gama baja-media y gama alta.

Para decidir el número óptimo de grupos que debemos crear, podemos usar la función NbClust de R, que nos devuelve cuál es (según unos criterios) el mejor número de clusters para el algoritmo de k-medias o bien, podemos ir probando con diferentes valores y decidir nosotros en función de la información que recabemos.

Primero usaremos la función, teniendo en cuenta que como máximo aceptaremos tener 10 grupos y como mínimo 2:

set.seed(220322)
nc = NbClust(TrainEscalado,min.nc=2,max.nc=10,method = "kmeans")

## *** : The Hubert index is a graphical method of determining the number of clusters.
##                 In the plot of Hubert index, we seek a significant knee that corresponds to a 
##                 significant increase of the value of the measure i.e the significant peak in Hubert
##                 index second differences plot. 
## 

## *** : The D index is a graphical method of determining the number of clusters. 
##                 In the plot of D index, we seek a significant knee (the significant peak in Dindex
##                 second differences plot) that corresponds to a significant increase of the value of
##                 the measure. 
##  
## ******************************************************************* 
## * Among all indices:                                                
## * 8 proposed 2 as the best number of clusters 
## * 3 proposed 3 as the best number of clusters 
## * 2 proposed 8 as the best number of clusters 
## * 11 proposed 9 as the best number of clusters 
## 
##                    ***** Conclusion *****                            
##  
## * According to the majority rule, the best number of clusters is  9 
##  
##  
## *******************************************************************
nc
## $All.index
##         KL       CH  Hartigan      CCC     Scott      Marriot   TrCovW
## 2   3.7372 2231.534  938.0726  -7.7989  5734.242 6.840925e+20 12493090
## 3   1.4122 1806.690  654.6000 -13.6059  9122.307 7.498264e+20  8370406
## 4   2.5713 1589.790  295.5897 -16.5777 13311.058 5.478838e+20  7359205
## 5   0.4411 1340.832  451.3764 -24.2826 15058.326 5.907891e+20  6585463
## 6  31.6152 1265.556  252.9865 -20.4498 16951.813 5.691667e+20  5176556
## 7   0.0315 1153.257  373.6215 -21.5790 18661.448 5.389226e+20  5009361
## 8   0.0906 1120.156 2916.5901 -17.5170 21486.427 3.864409e+20  4580569
## 9  93.9316 1952.131  139.3206  49.7614 32381.793 4.841428e+19  1449588
## 10  0.7282 1801.739  191.4079  46.5489 33128.923 5.100497e+19  1414613
##       TraceW Friedman  Rubin Cindex     DB Silhouette   Duda   Pseudot2   Beale
## 2  19173.790   8.5627 1.4739 0.0321 1.4470     0.4073 1.4789 -1208.2390 -1.2455
## 3  15988.705   9.7096 1.7675 0.0276 1.6252     0.3259 3.1466 -2281.9631 -2.6236
## 4  14037.001  13.5177 2.0133 0.0255 1.5512     0.3112 0.9610   100.4054  0.1562
## 5  13207.592  14.5017 2.1397 0.0249 1.6094     0.2273 1.6048  -815.8979 -1.4488
## 6  12051.656  15.4810 2.3449 0.0227 1.6606     0.2457 5.9409 -1353.1340 -3.1972
## 7  11436.708  17.0108 2.4710 0.0218 1.6345     0.2414 1.3219  -324.8093 -0.9358
## 8  10595.172  19.2072 2.6673 0.0210 1.5048     0.2557 2.3430  -772.6783 -2.1987
## 9   6539.603  29.1104 4.3214 0.0918 1.0827     0.2592 4.0564 -1253.0320 -2.8901
## 10  6351.410  29.9539 4.4494 0.0913 1.1398     0.2464 0.7969   190.8494  0.9793
##    Ratkowsky      Ball Ptbiserial    Frey McClain   Dunn Hubert SDindex Dindex
## 2     0.3323 9586.8949     0.4224  1.2264  0.4052 0.0012  1e-04  7.2473 1.5510
## 3     0.3600 5329.5683     0.4000  0.4116  0.8271 0.0016  1e-04  6.4723 1.3943
## 4     0.3326 3509.2502     0.4068  6.3164  0.9695 0.0007  1e-04  6.5693 1.2853
## 5     0.3074 2641.5184     0.3209  0.1855  1.7394 0.0011  1e-04  6.3678 1.2076
## 6     0.2932 2008.6094     0.3262  0.6156  1.9272 0.0003  1e-04  6.9963 1.1306
## 7     0.2772 1633.8154     0.3134 -0.0940  2.2187 0.0008  1e-04  6.6163 1.0874
## 8     0.2657 1324.3965     0.3223 -0.2051  2.1732 0.0008  1e-04  5.8757 1.0428
## 9     0.2920  726.6226     0.3274  4.3783  2.1327 0.0030  1e-04  6.6319 1.0283
## 10    0.2782  635.1410     0.3096  9.4740  2.4212 0.0030  1e-04  8.2962 1.0060
##      SDbw
## 2  1.3879
## 3  1.3910
## 4  1.1600
## 5  0.9978
## 6  1.0421
## 7  1.0252
## 8  0.8578
## 9  0.4463
## 10 0.7642
## 
## $All.CriticalValues
##    CritValue_Duda CritValue_PseudoT2 Fvalue_Beale
## 2          0.8638           588.2681       1.0000
## 3          0.8590           549.0712       1.0000
## 4          0.8575           410.6496       0.9879
## 5          0.8469           391.3096       1.0000
## 6          0.8456           297.1056       1.0000
## 7          0.8384           257.0461       1.0000
## 8          0.7997           337.5479       1.0000
## 9          0.7996           416.7944       1.0000
## 10         0.8391           143.6584       0.4373
## 
## $Best.nc
##                      KL       CH Hartigan     CCC    Scott      Marriot  TrCovW
## Number_clusters  9.0000    2.000    9.000  9.0000     9.00 9.000000e+00       3
## Value_Index     93.9316 2231.534 2777.269 49.7614 10895.37 3.406173e+20 4122684
##                   TraceW Friedman   Rubin Cindex     DB Silhouette   Duda
## Number_clusters    9.000   9.0000  9.0000  8.000 9.0000     2.0000 2.0000
## Value_Index     3867.376   9.9032 -1.5261  0.021 1.0827     0.4073 1.4789
##                  PseudoT2   Beale Ratkowsky     Ball PtBiserial   Frey McClain
## Number_clusters     2.000  2.0000      3.00    3.000     2.0000 2.0000  2.0000
## Value_Index     -1208.239 -1.2455      0.36 4257.327     0.4224 1.2264  0.4052
##                  Dunn Hubert SDindex Dindex   SDbw
## Number_clusters 9.000      0  8.0000      0 9.0000
## Value_Index     0.003      0  5.8757      0 0.4463
## 
## $Best.partition
##    [1] 4 4 9 7 1 5 9 9 6 7 4 9 1 5 7 9 6 6 4 1 7 6 5 9 7 1 9 5 1 9 9 1 9 9 9 2 5
##   [38] 4 7 7 9 4 6 2 1 7 2 6 6 6 4 4 9 6 4 5 6 9 7 4 1 1 9 7 9 7 5 4 5 7 7 9 5 4
##   [75] 6 7 5 6 4 6 6 6 5 6 9 6 9 4 6 4 6 9 7 9 4 5 9 9 6 5 5 6 4 9 4 6 7 5 9 5 5
##  [112] 4 4 7 5 2 4 6 5 6 9 9 6 4 6 9 7 9 2 6 9 9 4 7 9 7 1 4 5 4 1 9 4 6 5 5 9 9
##  [149] 6 1 7 9 6 7 7 9 9 4 6 7 4 1 6 7 4 1 4 9 6 5 4 9 9 7 4 4 1 4 1 9 5 6 1 5 6
##  [186] 4 2 4 3 5 6 4 6 6 9 1 9 2 5 6 7 9 7 1 6 6 7 1 4 6 6 6 6 4 5 9 9 9 5 5 9 6
##  [223] 1 9 9 4 2 6 4 6 9 9 3 6 5 9 1 7 4 9 7 7 6 4 5 6 7 9 5 7 9 7 4 6 7 5 5 9 2
##  [260] 6 7 5 6 6 6 5 6 9 2 5 6 2 9 4 9 1 6 6 5 4 9 5 9 6 1 5 2 6 7 4 1 4 6 4 9 7
##  [297] 5 9 6 1 9 6 6 1 9 5 6 4 6 6 1 9 5 4 6 1 9 9 5 7 5 2 4 7 5 6 5 1 9 6 7 7 9
##  [334] 9 7 9 6 9 3 5 1 9 4 5 5 4 6 5 6 1 9 7 6 5 1 7 5 9 5 9 9 4 9 6 9 7 6 6 5 1
##  [371] 4 4 1 7 3 5 9 9 5 7 2 7 6 1 7 5 7 5 5 6 6 4 5 5 5 7 9 5 6 9 9 9 6 2 7 7 2
##  [408] 6 4 9 6 2 7 7 7 9 7 7 7 9 5 5 7 7 4 9 5 5 7 6 2 5 9 5 9 6 5 1 9 9 9 6 9 6
##  [445] 4 5 7 6 6 5 6 9 7 6 5 5 6 2 7 1 9 6 9 6 5 7 4 4 9 4 4 5 6 7 5 6 4 6 6 4 1
##  [482] 5 9 2 5 4 6 4 2 1 9 4 7 7 5 5 7 2 6 7 7 9 4 6 7 7 5 9 9 9 9 4 5 3 6 6 4 5
##  [519] 7 7 6 9 5 9 4 5 9 7 7 6 9 7 5 6 6 9 6 2 5 5 4 4 4 7 7 5 9 5 7 1 6 2 9 6 4
##  [556] 5 5 6 9 7 9 6 6 5 6 7 9 9 4 1 1 5 6 9 4 2 6 5 1 5 9 4 2 5 7 6 5 5 7 6 6 5
##  [593] 9 5 6 5 6 5 9 3 2 9 5 4 9 5 5 6 4 6 7 5 9 6 4 6 4 6 6 9 4 4 6 6 9 6 2 7 6
##  [630] 9 9 9 7 4 6 9 6 6 3 1 7 9 4 9 9 5 5 6 7 5 1 7 4 6 7 9 7 6 2 1 4 6 1 9 7 4
##  [667] 2 5 2 6 9 5 2 9 6 4 4 5 5 6 1 9 5 5 5 6 5 5 4 7 6 1 9 4 4 5 6 5 6 7 5 9 1
##  [704] 6 9 4 6 6 5 6 9 6 4 6 5 9 2 5 9 6 4 9 6 6 7 5 3 7 5 9 7 5 6 5 1 5 9 5 9 1
##  [741] 9 6 4 7 8 4 4 7 9 7 9 6 3 9 9 6 5 6 5 1 7 1 4 9 5 5 9 6 5 4 6 5 7 9 9 6 6
##  [778] 1 6 9 6 6 9 9 5 5 2 9 9 9 6 5 4 6 9 4 9 2 4 6 7 6 9 5 5 5 7 5 5 1 5 7 9 9
##  [815] 9 2 7 6 5 5 4 6 6 2 9 1 9 6 5 9 6 6 9 4 5 6 6 6 5 7 5 9 5 4 6 9 9 6 5 9 5
##  [852] 9 9 1 5 5 4 9 6 6 9 6 7 6 6 4 7 1 9 6 4 7 5 9 6 7 6 7 1 9 9 5 5 5 6 7 7 5
##  [889] 6 4 7 9 4 4 6 2 1 1 5 7 5 5 9 5 9 4 9 1 7 7 6 5 6 6 7 7 7 6 5 4 4 9 5 7 9
##  [926] 4 6 5 4 4 9 7 9 6 9 5 6 5 6 9 1 6 5 7 5 7 9 6 7 5 9 6 5 9 7 6 5 7 7 6 6 6
##  [963] 9 7 4 4 5 5 6 4 5 9 9 6 5 5 9 5 4 5 4 6 7 9 7 4 9 6 5 7 7 5 7 3 6 7 4 6 9
## [1000] 7 2 5 7 4 5 5 4 9 5 9 7 1 7 7 1 9 9 6 6 9 9 9 6 9 5 4 6 2 7 5 9 6 6 5 5 6
## [1037] 6 5 7 9 5 1 9 9 7 5 6 7 2 5 6 6 4 9 6 4 9 6 2 9 1 9 6 7 6 9 6 6 5 1 3 9 9
## [1074] 7 7 6 7 5 1 5 9 6 6 9 6 4 9 7 5 6 6 7 6 6 6 1 6 1 3 7 9 4 2 5 7 9 9 6 6 3
## [1111] 9 7 5 9 9 7 9 9 4 5 9 6 9 1 6 5 9 9 6 1 7 4 6 7 7 1 5 6 6 7 5 6 6 1 1 9 9
## [1148] 5 5 6 9 6 7 5 6 4 1 7 9 6 9 9 6 6 7 7 7 9 9 9 9 6 6 1 1 6 6 9 9 9 5 5 3 4
## [1185] 5 7 9 9 3 2 1 5 5 5 9 4 7 9 1 9 6 2 4 5 4 9 6 9 1 5 5 1 9 9 4 5 4 6 4 9 2
## [1222] 6 7 1 9 7 5 6 6 5 6 2 4 5 2 6 1 5 9 6 4 6 4 4 3 6 5 9 7 6 7 1 4 6 6 6 6 7
## [1259] 7 5 9 5 6 2 1 9 9 9 1 5 9 9 5 3 1 4 6 9 7 3 6 1 9 6 4 6 6 6 5 9 4 6 4 4 7
## [1296] 1 2 7 9 7 9 5 9 5 5 1 5 5 6 9 5 5 9 5 9 5 5 5 1 6 9 1 6 1 7 9 4 5 9 6 6 5
## [1333] 4 5 5 7 5 7 4 1 1 7 5 6 9 5 9 9 5 6 6 6 7 9 5 6 6 6 5 4 5 4 9 5 7 6 7 4 9
## [1370] 4 4 9 7 6 6 6 7 7 9 1 6 9 5 2 9 1 1 9 7 9 9 7 5 9 4 4 1 9 1 7 1 4 4 6 7 5
## [1407] 4 5 7 7 5 7 6 9 7 1 9 7 6 6 7 9 9 9 5 4 1 6 9 6 7 4 9 9 6 6 3 9 9 5 7 4 9
## [1444] 5 6 2 5 7 9 7 9 7 7 9 6 7 9 4 6 6 9 2 5 5 6 5 5 9 2 1 9 1 5 7 5 4 7 9 4 2
## [1481] 1 9 5 1 5 9 6 6 9 9 7 6 3 4 5 3 6 5 6 4 6 5 5 3 4 6 7 6 9 2 9 9 6 6 9 4 7
## [1518] 5 5 6 5 7 5 9 4 9 7 7 6 9 6 7 6 6 7 7 9 5 6 7 6 6 5 1 5 1 7 9 7 9 9 1 7 6
## [1555] 1 6 6 5 5 5 7 1 5 9 1 9 7 7 1 6 9 5 9 5 5 2 2 6 7 4 9 6 5 6 3 5 6 9 9 9 7
## [1592] 6 4 5 2 6 6 9 5 4 9 1 7 6 9 4 4 4 6 1 9 5 6 6 2 1 7 6 9 4 9 9 9 9 9 9 7 1
## [1629] 5 2 4 9 1 5 6 6 9 6 4 9 5 9 5 9 9 5 6 5 5 6 6 9 5 6 6 5 9 4 2 9 6 5 5 1 6
## [1666] 4 5 4 7 6 5 4 1 9 7 6 4 6 9 7 5 6 7 7 9 9 9 1 5 6 9 2 9 6 7 9 7 6 7 6 6 6
## [1703] 6 7 4 6 9 5 6 9 6 2 5 2 9 4 4 9 6 6 7 1 4 1 5 7 6 5 5 6 6 4 5 4 9 7 9 5 5
## [1740] 7 4 6 6 7 9 7 5 9 4 2 6 3 9 5 7 7 5 7 5 5 2 9 9 4 7 5 4 5 9 6 6 9 4 1 9 9
## [1777] 5 6 9 7 5 4 3 9 9 5 1 4 6 6 9 9 4 5 1 9 4 2 4 6 9 7 7 9 7 6 5 2 9 5 9 6 7
## [1814] 9 5 9 5 5 5 6 7 9 9 1 4 4 5 4 1 4 5 1 7 1 7 6 5 7 7 5 6 6 4 9 7 6 4 6 4 5
## [1851] 4 9 9 9 6 1 9 5 4 7 1 4 9 5 9 1 6 9 4 5 6 6 9 1 6 9 9 5 9 4 9 6 5 6 7 9 2
## [1888] 1 7 6 6 7 1 6 9 5 5 9 4 4 7 9 6 9 5 9 6 9 6 6 6 4 4 3 4 6 7 4 6 1 7 7 4 6
## [1925] 9 6 6 6 7 7 9 1 2 4 1 5 2 9 6 3 1 6 4 6 7 4 6 6 7 6 5 2 4 9 6 3 4 4 4 5 6
## [1962] 6 2 1 4 7 9 1 7 7 4 4 7 6 9 9 5 7 9 5 4 4 5 4 4 1 6 5 5 9 6 6 9 7 7 4 7 7
## [1999] 9 6 4 9 5 4 4 6 6 6 5 9 5 5 7 9 7 9 6 4 6 6 1 5 9 1 9 6 9 5 6 9 6 4 4 9 4
## [2036] 7 6 9 9 9 6 7 6 6 9 9 9 1 1 7 6 4 6 9 4 5 6 9 9 5 7 4 6 9 9 6 6 4 4 6 4 7
## [2073] 4 7 9 4 3 5 6 4 4 4 4 7 6 9 9 5 2 5 5 5 9 5 4 9 6 9 9 5 4 7 4 6 6 5 7 9 9
## [2110] 6 3 6 5 9 6 5 6 9 5 5 6 5 4 4 6 6 7 7 9 9 6 7 9 5 9 9 5 4 7 6 5 5 5 6 2 6
## [2147] 1 5 9 7 7 5 7 6 9 2 5 7 4 4 1 4 1 6 9 6 9 9 9 5 4 9 9 6 6 5 2 6 4 2 5 7 9
## [2184] 1 9 3 1 4 5 2 1 4 6 1 7 9 7 5 6 5 9 1 7 5 9 4 4 1 6 7 4 5 1 6 6 7 1 5 9 6
## [2221] 4 7 7 9 7 5 9 9 5 9 9 6 4 6 7 6 7 7 7 6 6 6 9 4 4 6 6 5 6 6 6 4 5 1 5 9 9
## [2258] 5 9 5 6 7 9 5 9 9 7 5 9 1 4 6 9 1 4 6 9 6 7 6 7 9 5 7 7 7 5 4 1 7 7 6 9 7
## [2295] 6 2 4 6 5 5 1 5 6 7 4 9 6 4 9 6 6 7 9 9 4 5 6 9 7 4 9 4 7 9 9 7 9 6 9 5 6
## [2332] 2 5 9 5 4 3 6 9 6 9 3 7 6 7 5 5 4 1 6 5 6 4 6 6 6 9 9 5 4 4 5 3 9 4 5 6 9
## [2369] 5 5 6 4 9 4 5 5 5 5 4 9 4 2 7 7 6 5 6 4 5 4 1 5 7 6 6 4 7 7 6 9 9 9 6 1 5
## [2406] 6 9 6 9 6 6 5 1 7 9 7 9 9 3 3 6 9 6 5 4 5 9 6 5 6 1 4 1 6 4 9 9 7 9 9 9 7
## [2443] 6 9 9 7 7 5 6 5 6 4 5 9 6 1 9 5 1 5 7 9 1 9 7 7 1 4 7 1 4 4 5 5 6 4 4 9 7
## [2480] 6 2 6 6 5 6 6 7 3 7 4 7 9 5 6 9 4 9 2 7 7 9 7 5 7 4 9 7 5 6 6 9 6 6 5 5 5
## [2517] 5 4 6 6 6 6 3 7 5 6 5 9 6 5 1 7 5 9 7 4 9 7 4 7 6 4 7 4 1 9 9 5 7 9 9 6 4
## [2554] 4 4 5 6 9 5 4 6 5 9 6 5 9 7 7 9 4 4 5 4 7 9 2 6 6 2 1 6 9 7 5 9 7 4 5 5 6
## [2591] 4 6 5 5 7 4 9 4 6 6 9 9 7 1 5 7 9 6 5 5 5 7 7 5 2 7 7 9 7 7 9 9 7 9 6 7 6
## [2628] 7 1 1 9 4 5 9 9 7 5 5 9 9 9 6 6 5 5 6 4 6 4 5 1 5 4 5 9 6 7 5 6 2 6 6 6 6
## [2665] 9 6 6 7 9 9 1 5 6 7 9 4 5 9 6 9 7 4 5 6 4 6 6 5 6 5 5 9 5 6 9 6 6 6 5 6 9
## [2702] 1 7 2 1 4 4 6 4 4 9 5 4 5 9 7 6 7 5 1 5 2 9 6 1 4 1 6 6 9 6 6 6 1 6 4 9 6
## [2739] 1 9 4 5 2 9 5 9 9 4 9 3 7 9 7 5 7 6 5 5 5 7 9 5 6 9 6 1 6 2 4 6 6 5 7 6 1
## [2776] 4 1 4 2 5 7 9 4 6 4 6 5 2 9 4 9 6 6 5 6 2 2 9 9 6 9 3 7 9 7 4 7 5 6 5 2 5
## [2813] 9 9 6 2 6 9 6 9 5 9 4 6 9 5 6 1 6 4 9 5 6 5 1 7 6 9 6 7 5 7 7 9 2 9 4 1 9
## [2850] 9 4 7 7 7 5 9 7 5 1 5 9 5 9 5 6 5 7 6 7 5 9 7 4 9 4 6 7 5 7 9 6 4 7 6 7 6
## [2887] 9 6 5 7 4 7 5 9 7 9 9 9 6 7 3 6 6 9 5 6 9 4 7 9 9 5 9 6 4 9 5 6 5 4 5 5 3
## [2924] 4 6 7 7 9 6 9 1 9 9 9 5 7 4 6 1 4 6 6 9 7 2 4 9 7 5 2 5 4 5 5 6 4 4 7 9 9
## [2961] 5 1 1 9 7 5 5 5 4 7 6 6 5 7 1 5 1 6 6 7 5 6 6 1 4 5 4 9 6 4 9 9 4 6 9 9 7
## [2998] 2 7 9 9 9 9 7 2 5 9 5 9 5 9 4 5 7 7 6 1 5 5 1 7 5 6 5 4 1 2 4 5 7 9 7 6 9
## [3035] 4 7 7 4 4 1 5 5 4 7 4 6 7 9 4 7 7 6 7 9 1 6 6 7 6 5 7 5 7 1 5 7 9 1 6 1 4
## [3072] 6 3 9 2 7 9 9 6 5 9 5 1 9 7 9 6 4 7 6 9 7 5 9 7 5 5 6 2 6 5 5 7 5 9 1 7 4
## [3109] 5 7 6 6 4 7 4 7 7 9 9 6 5 7 7 9 6 4 5 5 9 1 6 7 9 7 6 7 9 5 9 9 9 4 1 4 9
## [3146] 1 6 6 6 7 1 5 6 6 9 6 1 6 7 9 9 9 9 5 4 5 6 7 7 1 6 9 4 4 4 7 9 1 6 6 6 5
## [3183] 9 7 4 9 4 1 7 6 5 5 5 7 6 9 4 9 9 9 5 9 9 6 6 4 6 3 7 7 7 9 5 1 4 9 4 9 6
## [3220] 1 5 4 7 7 9 9 5 6 6 4 4 6 7 5 5 9 5 9 4 6 1 3 1 4 4 2 2 5 9 5 9 4 7 6 1 4
## [3257] 1 6 6 9 9 7 6 3 4 7 7 4 4 9 7 5 5 4 4 9 7 9 1 5 4 6 3 5 3 6 6 9 9 4 4 4 5
## [3294] 5 4 5 9 1 7 6 7 2 6 1 5 9 5 9 4 9 5 1 5 9 1 7 1 7 9 7 5 7 7 9 6 4 4 6 5 5
## [3331] 6 7 4 7 7 4 5 9 7 9 7 9 7 6 9 9 1 7 9 5 5 5 9 7 6 1 6 4 7 6 6 5 7 4 6 6 9
## [3368] 5 1 5 9 5 4 5 9 1 5 7 9 6 6 9 4 7 9 2 7 6 4 2 1 7 7 7 7 2 9 9 9 5 5 1 4 4
## [3405] 4 1 9 5 9 6 6 5 5 7 1 6 7 9 4 9 5 7 9 5 6 6 6 7 2 6 5 9 9 1 6 7 2 9 9 7 6
## [3442] 5 7 6 7 5 5 7 4 7 4 4 7 7 3 4 9 1 4 5 7 9 9 6 4 4 5 6 9 4 6 9 6 4 9 4 2 6
## [3479] 9 6 7 5 6 9 4 9 3 6 4 7 4 4 9 5 6 9 6 5 4 6 6 6 4 5 6 6 9 6 7 9 5 5 9 7 4
## [3516] 6 7 9 5 7 7 1 7 4 2 6 5 1 6 7 9 5 9 6 6 9 5 2 7 9 6 4 6 6 2 7 5 6 9 9 6 4
## [3553] 6 5 9 5 6 2 5 4 9 7 6 9 9 6 9 6 6 9 6 7 6 5 6 5 6 9 9 7 4 7 5 6 6 4 2 6 6
## [3590] 3 5 5 4 6 5 1 7 7 5 9 4 4 2 2 4 4 7 9 4 9 5 5 7 4 2 9 6 4 6 7 1 6 4 6 7 5
## [3627] 4 5 4 7 6 4 1 9 1 9 5 9 1 6 7 6 9 6 5 3 7 5 2 7 5 9 5 7 7 7 4 5 2 6 5 9 6
## [3664] 9 9 9 1 7 5 5 2 5 5 4 5 4 6 9 7 7 7 5 2 9 1 4 6 5 9 6 6 4 1 7 9 5 9 7 9 7
## [3701] 9 1 5 6 7 9 9 7 9 9 4 6 3 6 6 4 5 5 6 7 7 6 9 5 7 6 9 5 1 5 7 5 7 2 9 9 6
## [3738] 7 4 9 5 6 7 4 9 5 9 9 9 5 2 1 9 4 9 6 2 2 6 7 1 5 2 6 7 9 4 9 6 9 7 7 6 9
## [3775] 9 9 5 5 7 9 7 5 4 4 5 6 6 9 9 5 2 5 6 4 7 6 5 5 6 4 1 6 6 9 5 6 2 5 7 7 5
## [3812] 4 9 4 6 4 5 6 4 5 7 4 9 9 6 4 1 4 6 5 1 5 9 1 9 6 5 7 5 4 4 9 2 5 9 6 4 6
## [3849] 6 6 1 5 5 5 2 6 9 5 6 4 9 9 7 1 7 1 5 5 2 2 4 9 5 6 1 6 1 1 4 6 9 9 9 5 9
## [3886] 9 4 9 6 6 5 4 1 9 6 6 9 9 5 9 5 1 9 6 6 7 9 7 7 4 6 9 5 4 4 9 2 5 9 2 9 9
## [3923] 4 9 6 4 9 6 4 7 9 5 7 9 6 3 4 4 5 6 2 3 5 4 6 1 4 9 1 9 4 5 4 7 2 2 9 4 7
## [3960] 5 9 1 6 9 7 1 4 1 5 5 9 5 6 4 7 6 7 9 9 6 6 9 6 9 6 1 4 6 5 5 4 2 2 7 4 6
## [3997] 9 5 7 6 7 7 5 9 6 1 5 5 9 9 4 6 9 9 5 6 4 6 5 9 6 6 5 6 9 5 6 9 6 5 9 5 9
## [4034] 5 4 6 9 1 9 5 7 3 1 5 4 9 4 4 1 1 6 4 1 5 5 5 1 6 6 7 7 9 6 5 9 5 2 5 5 5
## [4071] 9 5 6 6 2 4 6 6 6 6 6 5 9 6 9 9 6 7 4 4 9 4 1 1 7 5 6 6 4 7 6 9 9 9 6 6 5
## [4108] 9 5 7 9 2 9 2 7 5 9 2 6 1 6 5 5 7 4 5 9 7 7 9 5 6 3 5 9 9 5 6 9 9 9 5 5 6
## [4145] 9 7 5 7 6 5 7 6 7 5 6 2 4 6 9 9 9 9 9 9 7 6 6 5 7 5 9 9 2 7 1 6 9 5 1 5 4
## [4182] 4 9 9 9 6 9 5 5 9 4 9 7 9 1 9 1 6 6 2 1 9 9 4 4 4 9 7 6 4 6 6 4 7 5 5 6 6
## [4219] 5 9 6 5 6 5 6 9 4 9 1 6 7 9 4 7 7 6 4 5 7 9 9 7 9 5 4 6 7 5 6 1 5 4 9 7 7
## [4256] 7 7 6 5 7 7 7 4 9 6 7 6 9 6 4 5 4 7 6 6 6 4 1 5 2 7 5 1 9 5 1 7 9 6 6 9 6
## [4293] 6 6 7 2 7 3 4 9 7 5 5 6 6 6 4 6 9 9 6 5 1 6 9 2 9 6 4 2 6 5 4 5 1 5 4 7 4
## [4330] 7 1 1 4 4 4 9 7 7 6 7 7 6 6 9 6 7 5 4 6 6 7 5 9 4 9 9 5 7 4 5 4 3 5 6 6 4
## [4367] 4 4 4 9 6 9 4 5 6 4 6 1 1 6 1 9 9 6 4 4 5 6 4 1 5 6 1 6 5 9 6 7 5 9 6 5 6
## [4404] 7 1 6 4 5 6 1 9 7 5 7 9 9 5 5 9 9 9 6 6 9 7 6 9 1 5 5 5 5 9 6 9 7 3 6 1 6
## [4441] 9 6 5 5 9 5 5 5 9 7 6 9 2 6 4 1 6 1 7 2 6 4 9 5 9 7 6 5 1 6 4 7 2 1 4 1 9
## [4478] 5 5 9 9 7 4 6 6 7 5 5 7 7 9 6 5 6 4 6 6 6 7 1 6 6 9 9 5 6 5 6 5 5 7 5 9 5
## [4515] 5 9 9 5 9 7 1 9 9 9 1 6 7 7 4 9 4 9 7 5 1 9 6 3 5 7 5 7 9 9 9 5 7 7 9 7 6
## [4552] 9 6 5 9 9 9 1 6 7 7 9 5 4 5 6 9 9 9 9 4 6 6 3 6 7 1 6 7 1 4 4 7 9 1 6 5 4
## [4589] 7 3 6 7 6 9 5 3 6 2 1 6 7 2 9 6 6 5 6 5 6 7 1 4 5 5 5 6 6 6 4 7 6 2 1 5 9
## [4626] 6 6 9 4 6 9 2 9 9 4 7 6 6 9 9 1 9 9 7 5 4 6 6 3 6 6 6 4 4 9 2 5 6 2 5 4 3
## [4663] 7 1 4 9 5 7 7 6 1 5 9 9 6 9 9 9 6 6 4 9 6 6 9 9 7 4 9 9 6 4 9 9 5 9 4 5 5
## [4700] 2 4 7 4 9 7 6 9 6 6 5 9

La función NbClust prueba con diferente número de grupos y evalúa cuál es número óptimo de clusters según algunos criterios (muestra los resultados gráficos de la aplicación de algunos de ellos). Vemos que finalmente, nos dice que el número óptimo de grupos es 3, dado que es en el que más criterios de optimalidad coinciden.

Ahora, tras probar nosotros manualmente con diferente número de grupos, comprobamos que las mejores maneras para agrupar a los coches en función de la gama que tienen es creando 3,4 u 8 grupos.

# Ponemos una semilla para obtener siempre los mismos 3 grupos
set.seed(20112090)
km0 = kmeans(x=TrainEscalado, nstart = 5, centers = 3)
km0$centers #km2$[2]
##         Year Kilometers_Driven      Engine      Power      Price       kmpl
## 1  0.3253638        -0.0611996  1.37201471  1.4680576  1.4610785 -0.7413217
## 2 -1.0681538         0.2503815  0.04160836 -0.1492665 -0.4474343 -0.6201845
## 3  0.4455477        -0.1105182 -0.57602959 -0.5117040 -0.3478727  0.6340075
tkm0 <- table(km0$cluster,data$Gama)
tkm0
##    
##     Alta Baja-media
##   1  905         73
##   2  356        953
##   3   92       2332
# Ponemos una semilla para obtener siempre los mismos 4 grupos
set.seed(22032023)

km1 = kmeans(x=TrainEscalado, nstart = 5, centers = 4)
km1$centers #km1$[2]
##         Year Kilometers_Driven     Engine      Power      Price       kmpl
## 1 -1.2013643         0.1464077 -0.4491522 -0.5128053 -0.5978962 -0.3443545
## 2 -0.1286413         0.1717612  1.0387989  0.8212333  0.3662194 -0.7705951
## 3  0.5202489        -0.1219503 -0.5769438 -0.5043002 -0.3395386  0.6622250
## 4  0.6284743        -0.1958552  1.7519820  2.1655859  2.8047656 -0.7944150
tkm1 <- table(km1$cluster,data$Gama)
tkm1
##    
##     Alta Baja-media
##   1   77        953
##   2  870        232
##   3   74       2159
##   4  332         14
# Ponemos una semilla para obtener siempre los mismos 8 grupos
set.seed(20112020)
km2 = kmeans(x=TrainEscalado, nstart = 5, centers = 8)
km2$centers #km2$[2]
##         Year Kilometers_Driven     Engine      Power       Price       kmpl
## 1  0.7474046      -0.215013813 -0.4558985 -0.3368124 -0.23872151  0.1001606
## 2  0.6705904      -0.143829782  0.9138645  1.0616096  1.31598716 -0.4591534
## 3  0.6179466      -0.089033566 -0.7068906 -0.6775746 -0.38607087  1.3429620
## 4 -0.1313711       0.006899268 -0.7771851 -0.8573284 -0.47245641 -3.8415039
## 5 -0.5355855       0.062374606 -0.5472943 -0.5250797 -0.53531404  0.2572020
## 6 -2.1245518       0.281847224 -0.3699003 -0.5410684 -0.67361109 -0.3691450
## 7  0.2505949      -0.142263332  2.5613160  3.0441732  3.36430324 -1.0824906
## 8 -0.7006775       0.367237420  1.1640519  0.7542959 -0.01418397 -0.9321174
table(km2$cluster,data$Gama)
##    
##     Alta Baja-media
##   1   50        977
##   2  545         64
##   3   22        803
##   4    5         61
##   5   29        977
##   6   32        300
##   7  166          2
##   8  504        174
# Ponemos una semilla para obtener siempre los mismos 8 grupos
set.seed(201120572)
km3 = kmeans(x=TrainEscalado, nstart = 5, centers = 9)
km3$centers #km2$[2]
##         Year Kilometers_Driven     Engine      Power       Price        kmpl
## 1  0.7457212      -0.216085906 -0.4606091 -0.3420833 -0.24078272  0.08982149
## 2 -0.5380988       0.062578193 -0.5479289 -0.5259096 -0.53567271  0.25766538
## 3  0.2505949      -0.142263332  2.5613160  3.0441732  3.36430324 -1.08249057
## 4  0.6744025      -0.143545245  0.9136031  1.0614304  1.32163596 -0.45638120
## 5 -0.6978533       0.270519256  1.1603176  0.7529817 -0.01917599 -0.93074958
## 6 -0.1313711       0.006899268 -0.7771851 -0.8573284 -0.47245641 -3.84150392
## 7  1.1117978      63.738691629  2.2760090  2.7023168  4.98013063 -0.42798303
## 8  0.6197888      -0.088702794 -0.6984679 -0.6671467 -0.38233324  1.33768775
## 9 -2.1245518       0.281847224 -0.3699003 -0.5410684 -0.67361109 -0.36914499
table(km3$cluster,data$Gama)
##    
##     Alta Baja-media
##   1   50        968
##   2   29        973
##   3  166          2
##   4  542         63
##   5  506        176
##   6    5         61
##   7    1          0
##   8   22        815
##   9   32        300

Para no complicar demasiado el entendimiento del algoritmo, decidimos quedarnos con 3 clusters:

  • En el primer cluster, sobre todo, clasifica a los coches con valores altos de Price, Engine y Power. En esta categoría tenemos 978 coches, es decir, al 20.76% del total.

  • En el segundo cluster, si nos fijamos, clasifica a los coches con valores medios en casi todas las variables, con año de fabricación antiguo y kmpl alto. En esta categoría tenemos 1309 coches, es decir, al 27.79% del total.

  • En el tercer cluster, clasificamos a los coches con valores bajos en Precio, Engine y Power. En esta categoría tenemos 2424 coches relativamente comunes (sin valores atípicos o muy influyentes en ninguna de sus variables), es decir, al 51.45% del total.

Fijándonos en la tabla que nos devuelve, vemos que en el primer grupo la mayoría de los coches son de gama alta (el 92.54%); en el segundo grupo, los coches están bastante mezclados aunque son mayoría en la clase de gama media-baja pero necesitaríamos analizarlos más en profundidad para poder separarlos mejor (72.8% de coches de gama alta frente a un 27.2% de coches de gama media o baja); y en el tercer grupo, al contario que en el primero, la mayoría son de gama baja o media (un 96.2%)

Veamos esto que acabamos de explicar con algunos gráficos.

Según los cluster que hemos formado, si enfrentamos el precio (dominante en el grupo 1 con valores altos) frente a la variable Engine (dominante en el grupo 3 con valores bajos), deberíamos obtener un gráfico en el que los cluster 1 y 3 estuviesen bien diferenciados y el 2, que tenía valores medios, esté mezclado con ambos.

plot(TrainEscalado$Engine, TrainEscalado$Price, col=km0$cluster, pch=19 , cex=2, xlab = "Engine", ylab="Price", main = "Engine vs Price")
legend(-2,12,c("Cluster 1","Cluster 2","Cluster 3"),fill = (unique(sort(km0$cluster))))

Observamos como los coches de los cluster 1 y 3 están perfectamente separados. Además, esperábamos que el cluster 3 estuviese mezclado con los otros dos, sin embargo, vemos que enfrentando estas dos variables, separamos muy bien los tres grupos pese a que si que se juntan en ciertos coches.

Podríamos seguir haciendo gráficos y comprobaciones, pero con eso ya vemos que tenemos una buena forma de clasificar a algunos de los coches. En el grupo 1 teníamos mayoritariamente coches de gama alta; en el grupo 2, mezcla; y en el grupo 3, coches de gama baja o media. Hemos comprobado que los grupos 1 y 3 se separan muy bien gráficamente, pero de hecho, los grupos 1 y 2 también, dado que donde más "mezcla apreciamos es entre los cluster 2 y 3. Esto nos es de gran ayuda.

Clustering jerárquico

Para realizar el clustering jerárquico, también utilizaremos la distancia euclídea. Pero debemos definir cuál será la distancia entre dos grupos, que será la que nos sirva como criterio para decidir cuando se deben unir dos cluster. En R podemos definir diferentes distancias entre ellos (distancias entre los centroides de cada grupo, distancias entre los elementos más próximos de cada grupo, distancias entre los elementos más alejados de cada grupo…). Veremos cuál es el resultado de utilizar alguna de ellas mostrando los dendrogramas asociados a cada una.

En el primer caso utilizaremos el método single (distancia entre los elementos más cercanos de cada cluster):

#TrainEscalado[-745,] <- TrainEscalado
hc1 = hclust(d=dist(TrainEscalado), method = "single" )
plot(hc1)

En este dendrograma, podemos ver que tenemos un coche atípico (el 745), que se acumula a la izquierda.

En el segundo caso utilizaremos el método complete (distancia entre los elementos más alejados de cada cluster):

hc2 = hclust(d=dist(TrainEscalado), method = "complete")
plot(hc2)

En este caso sucede algo parecido a lo que ocurría antes, el coche 745, vuelve a aparecer “al margen” del resto, y vuelve a destacar por ser uno de los últimos coches en agruparse en el dendrograma, aunque ya se empiezan a visualizar diferentes agrupaciones.

En el tercer caso utilizaremos el método average (distancia media entre las observaciones de cada cluster):

hc3 = hclust(d=dist(TrainEscalado), method = "average" )
plot(hc3)

En este tercer caso, de nuevo, volvemos a observar que el coche 745 está más separado del resto que los demás entre sí, y que es este coche el que provoca el “retraso” de la unión de las diferentes agrupaciones iniciales.

En el cuarto y último caso utilizaremos el método centroid (distancia entre los centroides de cada cluster):

hc4 = hclust(d=dist(TrainEscalado), method = "centroid" )
plot(hc4)

De nuevo, volvemos a identificar al coche 745 a la izquierda, y volvemos a comprobar cómo es este coche el que produce un mayor retraso en la última unión de todos los grupos.

Como ya se vio en el análisis exploratorio de datos, podríamos pensar en eliminarlo, pero como tampoco tenemos más información y no sabemos si en realidad son outliers o simplemente tienen características un poco diferentes a las de los demás, decidimos quedarnos con ellos, ya que puede que nos aparezca más adelante otro coche similar a alguno de ellos, y nos ayuden a clasificarlo adecuadamente.

Por último, tomando el segundo dendrograma, que parece ser en el que mejor se observan las diferentes agrupaciones, vamos a proceder a tomar únicamente 3 grupos (igual que hicimos en el clustering no jerárquico).

plot(hc2)
rect.hclust(hc2, k=4, border = "blue")

Con esto vemos que tenemos un cluster principal en el que se aglutinan la mayoría de los coches y otros dos formados por muy pocos coches (los más atípicos de los extremos del dendrograma). Veamos cómo al formar 10 grupos se observan mejor la subdivisiones.

plot(hc2)
rect.hclust(hc2, k=10, border = "blue")

Análisis supervisado

Regresión logística

La división de los coches en gama baja-media y gama alta ha sido necesaria también para poder realizar la regresión logística sobre la nueva variable.

Se realizan las siguientes gráficas para visualizar los datos con la nueva variable:

table(data$Gama)
## 
##       Alta Baja-media 
##       1353       3358
ggplot(data, aes(x = Power, y = Price, color = Gama)) + geom_point()

ggplot(data, aes(x = Engine, y = Price, color = Gama)) + geom_point()

ggplot(data, aes(x = Transmission, y = Price, color = Gama)) + geom_point()

ggplot(data, aes(x = Seats, y = Price, color = Gama)) + geom_point()

ggplot(data, aes(x = Owner_Type, y = Price, color = Gama)) + geom_point()

Ahora, se van a aplicar los modelos con distintas covariables para buscar el mejor de ellos:

# modelos lineales generalizados estimados por MLE
logit <- glm( 
  Gama ~Power+Seats+Transmission+Owner_Type+Engine, 
  data = data, 
  family = binomial()
)
summary(logit)
## 
## Call:
## glm(formula = Gama ~ Power + Seats + Transmission + Owner_Type + 
##     Engine, family = binomial(), data = data)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.9207  -0.0737   0.1778   0.3066   4.0002  
## 
## Coefficients:
##                            Estimate Std. Error z value Pr(>|z|)    
## (Intercept)              -3.668e+00  4.878e+02  -0.008   0.9940    
## Power                    -1.940e-02  2.370e-03  -8.183 2.76e-16 ***
## Seats4                    1.186e+01  4.878e+02   0.024   0.9806    
## Seats5                    1.236e+01  4.878e+02   0.025   0.9798    
## Seats6                   -6.639e+00  6.166e+02  -0.011   0.9914    
## Seats7                    1.158e+01  4.878e+02   0.024   0.9811    
## Seats8                    9.311e+00  4.878e+02   0.019   0.9848    
## Seats9                    1.248e+01  4.878e+02   0.026   0.9796    
## Seats10                   1.171e+01  4.878e+02   0.024   0.9809    
## TransmissionManual        8.966e-01  1.476e-01   6.073 1.25e-09 ***
## Owner_TypeSecond         -6.252e-02  1.493e-01  -0.419   0.6755    
## Owner_TypeThird           6.359e-01  3.765e-01   1.689   0.0912 .  
## Owner_TypeFourth & Above  1.862e+00  1.113e+00   1.672   0.0945 .  
## Engine                   -3.224e-03  2.365e-04 -13.630  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 5649.7  on 4710  degrees of freedom
## Residual deviance: 2164.7  on 4697  degrees of freedom
## AIC: 2192.7
## 
## Number of Fisher Scoring iterations: 15
logit2 <- glm( 
  Gama ~Owner_Type+Seats+Transmission+Engine, 
  data = data, 
  family = binomial()
)
summary(logit2)
## 
## Call:
## glm(formula = Gama ~ Owner_Type + Seats + Transmission + Engine, 
##     family = binomial(), data = data)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.8689  -0.0929   0.1818   0.3542   4.0487  
## 
## Coefficients:
##                            Estimate Std. Error z value Pr(>|z|)    
## (Intercept)              -5.076e+00  4.895e+02  -0.010   0.9917    
## Owner_TypeSecond         -1.027e-02  1.475e-01  -0.070   0.9445    
## Owner_TypeThird           7.748e-01  3.781e-01   2.049   0.0405 *  
## Owner_TypeFourth & Above  2.289e+00  1.152e+00   1.987   0.0469 *  
## Seats4                    1.256e+01  4.895e+02   0.026   0.9795    
## Seats5                    1.327e+01  4.895e+02   0.027   0.9784    
## Seats6                   -5.571e+00  6.100e+02  -0.009   0.9927    
## Seats7                    1.297e+01  4.895e+02   0.026   0.9789    
## Seats8                    1.077e+01  4.895e+02   0.022   0.9824    
## Seats9                    1.508e+01  4.895e+02   0.031   0.9754    
## Seats10                   1.423e+01  4.895e+02   0.029   0.9768    
## TransmissionManual        1.361e+00  1.326e-01  10.269   <2e-16 ***
## Engine                   -4.555e-03  1.894e-04 -24.056   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 5649.7  on 4710  degrees of freedom
## Residual deviance: 2231.2  on 4698  degrees of freedom
## AIC: 2257.2
## 
## Number of Fisher Scoring iterations: 15
logit3 <- glm( 
  Gama ~Power+Transmission+Engine, 
  data = data, 
  family = binomial()
)
summary(logit3)
## 
## Call:
## glm(formula = Gama ~ Power + Transmission + Engine, family = binomial(), 
##     data = data)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.8679  -0.0698   0.1875   0.3241   4.2115  
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         9.3942879  0.3385829  27.746  < 2e-16 ***
## Power              -0.0140166  0.0021605  -6.488 8.73e-11 ***
## TransmissionManual  0.6373511  0.1373724   4.640 3.49e-06 ***
## Engine             -0.0040538  0.0001737 -23.345  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 5649.7  on 4710  degrees of freedom
## Residual deviance: 2298.2  on 4707  degrees of freedom
## AIC: 2306.2
## 
## Number of Fisher Scoring iterations: 6

La interpretación de los p-valores es similar a la del modelo lineal. Podemos ver que las variables Engine,Power y Transmission son significativas en el modelo (p-valor mucho menor de 0.05), mientras que la variable Seats y Owner_Type influyen más en un modelo que en otro.

El mejor modelo es el explicado por las variables Power, Engine y Transmission. En cuanto a los coeficientes, la interpretación cambia. El modelo GLM no ajusta la variable respuesta sino una función de enlace. En el caso del modelo logit esta función es: \(η=log(p1−p)\), siendo \(p\) la probabilidad de que el individuo tome el valor “1” correspondiente a la gama alta en la variable dicotómica. Al cociente \(p/(1−p)\) se le conoce como odds ratio. Por tanto, los coeficientes del modelo logit se interpretan como el logaritmo del odds ratio. Si nos fijamos en el coeficiente de la variable Transmission (0.63) en el modelo 3, nos está indicando que el logaritmo del odds ratio de pertenecer al grupo de los coches de alta gama aumenta 0.63 unidades por cada unidad que aumenta la variable Transmission.

Antes de comenzar con las siguientes, lo que debemos hacer es definir una medida de precisión para contrastar los datos una vez que tengamos cada matriz de confusión y comparar los resultados que nos ofrecen cada uno de los métodos que empleemos. En nuestro caso, la variable respuesta no está muy bien balanceada:

summary(data$Gama)
##       Alta Baja-media 
##       1353       3358

Como vemos, los coches de gama alta se corresponden con el 30% aproximadamente y los de gama baja-media, el 70%. Supondremos que están más o menos equilibradas. No queremos dar más valor a identificar un tipo de coche frente al otro. Por tanto, se establecerá la siguiente medida de precisión. Esta es una medida de precisión que hemos creado para tener en cuenta todos los casos, tanto los falsos negativos como los falsos positivos. En realidad, se trata de la media geométrica de la sensitividad (recall) y la especificidad, y se define según la siguiente expresión:

\[\text{Medida de precisión}=\sqrt{\frac{TP}{FN+TP}·\frac{TN}{FP+TN}}=\sqrt{TPR·TNR} \] donde:

  • \(TP\) (true positive) son los coches de gama alta que acertamos que son de gama alta.

  • \(TN\) (true negative) son los coches de gama baja-media que acertamos que son de baja-media.

  • \(FP\) (false positive) son los coches de gama baja-media que nosotros predecimos como gama alta.

  • \(FN\) (false negative) son los coches de gana alta que nosotros predecimos como de gama baja-media.

  • \(TPR\) (sensitivity, recall, hit rate or true positive rate) es la sensitividad.

  • \(TNR\) (specificity, dplyr::selectivity or true negative rate) es la especificidad.

Cabe señalar que esta medida sólo será utilizada en futuros análisis cuando el problema requiera de una precisión de este tipo, mientras, se utilizará como métrica el accuracy.

medidaPrecision <- function(matrizDeConfusion){
  return(sqrt(matrizDeConfusion[1,1]*matrizDeConfusion[2,2]/((matrizDeConfusion[1,2]+matrizDeConfusion[2,2])*(matrizDeConfusion[2,1]+matrizDeConfusion[1,1]))))
}

Los k-vecinos más cercanos (k-NN: The k-nearest neighbours)

El método de los k-vecinos más cercanos consiste en clasificar a un nuevo individuo en función de la categoría de sus \(k\) vecinos más cercanos, es decir, clasificaremos un coche en gama alta o gama media-baja en función de la gama de los coches más cercanos a él (con cercanía nos referimos a similitud entre sus características).

# Ponemos una semilla para que siempre nos salga el mismo resultado (el algoritmo es aleatorio)
set.seed(07122020)
trainX <- data[,c(2:10,13,14)]
preProcValues <- preProcess(x = trainX,method = c("center", "scale"))
preProcValues
## Created from 4711 samples and 11 variables
## 
## Pre-processing:
##   - centered (6)
##   - ignored (5)
##   - scaled (6)

Para poder realizar las predicciones, se entrena mediante cross-validation y así se estima el número óptimo de vecinos.

set.seed(400)
ctrl <- trainControl(method="repeatedcv",repeats = 3)
knnFit <- train(Gama ~ ., data = data[,c(2:10,13,14,15)], method = "knn", trControl = ctrl, preProcess = c("center","scale"), tuneLength = 20)
knnFit
## k-Nearest Neighbors 
## 
## 4711 samples
##   11 predictor
##    2 classes: 'Alta', 'Baja-media' 
## 
## Pre-processing: centered (30), scaled (30) 
## Resampling: Cross-Validated (10 fold, repeated 3 times) 
## Summary of sample sizes: 4240, 4240, 4240, 4240, 4240, 4240, ... 
## Resampling results across tuning parameters:
## 
##   k   Accuracy   Kappa    
##    5  0.9133979  0.7882405
##    7  0.9133945  0.7892998
##    9  0.9110582  0.7844472
##   11  0.9080157  0.7778724
##   13  0.9053970  0.7721609
##   15  0.9045491  0.7700798
##   17  0.9045492  0.7699856
##   19  0.9049044  0.7700120
##   21  0.9041972  0.7677520
##   23  0.9049056  0.7685929
##   25  0.9047639  0.7677360
##   27  0.9049062  0.7678410
##   29  0.9044819  0.7665081
##   31  0.9039856  0.7648794
##   33  0.9029251  0.7617441
##   35  0.9024993  0.7604944
##   37  0.9015801  0.7579416
##   39  0.9023573  0.7594855
##   41  0.9024273  0.7591610
##   43  0.9029213  0.7599820
## 
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was k = 5.
plot(knnFit)

knnPredict <- predict(knnFit,newdata =dataTest[,c(2:10,13,14,15)])
data$Gama <- factor(data$Gama,levels=c("Alta","Baja-media"),labels=c("Alta","Baja-media"))
dataTest$Gama <- factor(dataTest$Gama,levels=c("Alta","Baja-media"),labels=c("Alta","Baja-media"))
confusionMatrix(knnPredict, dataTest$Gama )
## Confusion Matrix and Statistics
## 
##             Reference
## Prediction   Alta Baja-media
##   Alta        275         39
##   Baja-media   61        798
##                                           
##                Accuracy : 0.9147          
##                  95% CI : (0.8973, 0.9301)
##     No Information Rate : 0.7136          
##     P-Value [Acc > NIR] : < 2e-16         
##                                           
##                   Kappa : 0.7873          
##                                           
##  Mcnemar's Test P-Value : 0.03573         
##                                           
##             Sensitivity : 0.8185          
##             Specificity : 0.9534          
##          Pos Pred Value : 0.8758          
##          Neg Pred Value : 0.9290          
##              Prevalence : 0.2864          
##          Detection Rate : 0.2344          
##    Detection Prevalence : 0.2677          
##       Balanced Accuracy : 0.8859          
##                                           
##        'Positive' Class : Alta            
## 
accuracy <- mean(knnPredict == dataTest$Gama)
accuracy
## [1] 0.9147485
error = 1-accuracy
error
## [1] 0.08525149
table(knnPredict,dataTest$Gama)
##             
## knnPredict   Alta Baja-media
##   Alta        275         39
##   Baja-media   61        798

Por tanto, tomando el número óptimo de vecinos para realizar las predicciones, tenemos una precisión del 91’5%.

Árboles de decisión

Los árboles de decisión se construyen en base al cumplimiento o no de ciertos criterios en torno a las variables explicativas de los datos. Digamos que comenzamos comprobando un cierto criterio sobre la variable más explicativa. Si el criterio se cumple, nos vamos por una rama; si no, nos vamos por la otra. Por cada una de ellas, volveremos a comprobar otro cierto criterio, y así hasta llegar a las hojas. Las hojas clasifican a los datos en un grupo u otro. Todo esto lo explicaremos mucho mejor cuando construyamos nuestros árboles de decisión.

Primero veremos cuál es el valor óptimo de cp (parámetro de complejidad que combina la tasa de error con la profundidad del árbol). Para ello, construiremos un árbol muy complejo, y en función de los resultados que obtengamos lo “podaremos”.

arbolDecisionInfoComplejo <- rpart(data[,15] ~., data[,c(2:10,13,14)], cp=0.001,  method="class",parms=list(split="information"))
plotcp(arbolDecisionInfoComplejo)

En el gráfico vemos que el mínimo se encuentra en torno al 0,0036, así que para obtener el mejor compromiso entre error y profundidad nos quedaremos con un valor de cp=0,004, porque tenemos que tener en cuenta que cuanto más pequeño sea este valor, más complejo y profundo será el árbol y menos error tendremos, sin embargo, corremos el riesgo de sobreajustar el modelo (problema de overfitting). Por otro lado, podemos usar dos criterios diferentes: el de información y el de Gini. Usaremos los dos y comprobaremos con cuál obtenemos una mayor precisión.

Comenzamos con el árbol de decisión utilizando el criterio de información.

arbolDecision <- rpart(data[,15] ~., data[,c(2:10,13,14)], cp=0.004,  method="class",parms=list(split="information"))
arbolDecision
## n= 4711 
## 
## node), split, n, loss, yval, (yprob)
##       * denotes terminal node
## 
##   1) root 4711 1353 Baja-media (0.287200170 0.712799830)  
##     2) Engine>=1690 1548  291 Alta (0.812015504 0.187984496)  
##       4) Power>=164.2 860   65 Alta (0.924418605 0.075581395)  
##         8) Fuel_Type=Diesel 726   31 Alta (0.957300275 0.042699725)  
##          16) Seats=4,5,6,8 533    7 Alta (0.986866792 0.013133208) *
##          17) Seats=7 193   24 Alta (0.875647668 0.124352332)  
##            34) Engine>=2299.5 177   10 Alta (0.943502825 0.056497175)  
##              68) kmpl>=10.955 166    0 Alta (1.000000000 0.000000000) *
##              69) kmpl< 10.955 11    1 Baja-media (0.090909091 0.909090909) *
##            35) Engine< 2299.5 16    2 Baja-media (0.125000000 0.875000000) *
##         9) Fuel_Type=Petrol 134   34 Alta (0.746268657 0.253731343)  
##          18) Power>=180.515 96    9 Alta (0.906250000 0.093750000) *
##          19) Power< 180.515 38   13 Baja-media (0.342105263 0.657894737)  
##            38) Power< 177.235 7    0 Alta (1.000000000 0.000000000) *
##            39) Power>=177.235 31    6 Baja-media (0.193548387 0.806451613) *
##       5) Power< 164.2 688  226 Alta (0.671511628 0.328488372)  
##        10) Engine>=2117.5 420   68 Alta (0.838095238 0.161904762)  
##          20) Power< 140.5 335   22 Alta (0.934328358 0.065671642)  
##            40) Power< 137 229    6 Alta (0.973799127 0.026200873) *
##            41) Power>=137 106   16 Alta (0.849056604 0.150943396)  
##              82) Power>=138.555 90    0 Alta (1.000000000 0.000000000) *
##              83) Power< 138.555 16    0 Baja-media (0.000000000 1.000000000) *
##          21) Power>=140.5 85   39 Baja-media (0.458823529 0.541176471)  
##            42) kmpl>=13.59 31    6 Alta (0.806451613 0.193548387) *
##            43) kmpl< 13.59 54   14 Baja-media (0.259259259 0.740740741)  
##              86) Engine< 2188.5 9    0 Alta (1.000000000 0.000000000) *
##              87) Engine>=2188.5 45    5 Baja-media (0.111111111 0.888888889) *
##        11) Engine< 2117.5 268  110 Baja-media (0.410447761 0.589552239)  
##          22) Engine< 1796.5 22    0 Alta (1.000000000 0.000000000) *
##          23) Engine>=1796.5 246   88 Baja-media (0.357723577 0.642276423)  
##            46) Power< 155.875 210   88 Baja-media (0.419047619 0.580952381)  
##              92) Engine< 1798.5 28    1 Alta (0.964285714 0.035714286) *
##              93) Engine>=1798.5 182   61 Baja-media (0.335164835 0.664835165)  
##               186) Engine>=1932 134   61 Baja-media (0.455223881 0.544776119)  
##                 372) kmpl< 16.88 76   26 Alta (0.657894737 0.342105263)  
##                   744) kmpl>=15.05 39    2 Alta (0.948717949 0.051282051) *
##                   745) kmpl< 15.05 37   13 Baja-media (0.351351351 0.648648649) *
##                 373) kmpl>=16.88 58   11 Baja-media (0.189655172 0.810344828)  
##                   746) kmpl>=20.19 8    0 Alta (1.000000000 0.000000000) *
##                   747) kmpl< 20.19 50    3 Baja-media (0.060000000 0.940000000) *
##               187) Engine< 1932 48    0 Baja-media (0.000000000 1.000000000) *
##            47) Power>=155.875 36    0 Baja-media (0.000000000 1.000000000) *
##     3) Engine< 1690 3163   96 Baja-media (0.030350933 0.969649067)  
##       6) Engine>=1353.5 1250   80 Baja-media (0.064000000 0.936000000)  
##        12) Engine< 1366 30    0 Alta (1.000000000 0.000000000) *
##        13) Engine>=1366 1220   50 Baja-media (0.040983607 0.959016393)  
##          26) Seats=4,7 59   19 Baja-media (0.322033898 0.677966102)  
##            52) Power>=99.41 31   13 Alta (0.580645161 0.419354839)  
##             104) Fuel_Type=Diesel 15    0 Alta (1.000000000 0.000000000) *
##             105) Fuel_Type=Petrol 16    3 Baja-media (0.187500000 0.812500000) *
##            53) Power< 99.41 28    1 Baja-media (0.035714286 0.964285714) *
##          27) Seats=5,8 1161   31 Baja-media (0.026701120 0.973298880)  
##            54) Engine< 1496.5 416   27 Baja-media (0.064903846 0.935096154)  
##             108) Engine>=1495.5 18    0 Alta (1.000000000 0.000000000) *
##             109) Engine< 1495.5 398    9 Baja-media (0.022613065 0.977386935) *
##            55) Engine>=1496.5 745    4 Baja-media (0.005369128 0.994630872) *
##       7) Engine< 1353.5 1913   16 Baja-media (0.008363826 0.991636174)  
##        14) Seats=6 7    0 Alta (1.000000000 0.000000000) *
##        15) Seats=4,5,7,8 1906    9 Baja-media (0.004721931 0.995278069) *

Aquí tenemos nuestro árbol de decisión. Como vemos, la primera pregunta (2 y 3) que le hará a un nuevo coche es si su motor es mayor o igual o menor a 1690. Si es mayor o igual a 1690, le preguntará si la potencia es inferior a 164.2 (preguntas 4 y 5); y si es que no, preguntará si el tipo de combustible es de diésel o petróleo(preg 8 y 9). Si es de diésel, entonces habrá que fijarse en el número de asientos que posea el coche. En este punto del árbol, si el número de asientos es 7, se volverá a clasificar por el motor, en función de si es mayor a 2299 en el que cerraremos la rama del árbol observando si la variable kmpl es inferior a 10. Así se ha llegado a un nodo hoja (marcados en asterisco) como por ejemplo llegar hasta la condición 69 y cumplirla, donde se clasifica el coche en gama media-baja o gama alta (dependiendo de la hoja en cuestión a la que haya llegado el coche). Como acabamos de señalar, las variables más informativas de nuestro árbol son el motor y la potencia de los coches. Tiene sentido que esto sea así, y que la primera característica sobre la que “pregunte” el árbol para clasificar a un nuevo coche sea Engine. ¿Por qué? Porque tal y como vimos, en la sección del análisis exploratorio de datos, justo era ésta la covariable más relacionada con gama de los coches. Así que los resultados obtenidos son perfectamente coherentes con todo el estudio realizado hasta el momento.

Aquí tenemos las preguntas que se han hecho a lo largo del árbol.

labels(arbolDecision, pretty=T)
##  [1] "root"           "Engine>=1690"   "Power>=164.2"   "Fuel_Type=Disl"
##  [5] "Seats=4,5,6,8"  "Seats=7"        "Engine>=2300"   "kmpl>=10.96"   
##  [9] "kmpl< 10.96"    "Engine< 2300"   "Fuel_Type=Ptrl" "Power>=180.5"  
## [13] "Power< 180.5"   "Power< 177.2"   "Power>=177.2"   "Power< 164.2"  
## [17] "Engine>=2118"   "Power< 140.5"   "Power< 137"     "Power>=137"    
## [21] "Power>=138.6"   "Power< 138.6"   "Power>=140.5"   "kmpl>=13.59"   
## [25] "kmpl< 13.59"    "Engine< 2188"   "Engine>=2188"   "Engine< 2118"  
## [29] "Engine< 1796"   "Engine>=1796"   "Power< 155.9"   "Engine< 1798"  
## [33] "Engine>=1798"   "Engine>=1932"   "kmpl< 16.88"    "kmpl>=15.05"   
## [37] "kmpl< 15.05"    "kmpl>=16.88"    "kmpl>=20.19"    "kmpl< 20.19"   
## [41] "Engine< 1932"   "Power>=155.9"   "Engine< 1690"   "Engine>=1354"  
## [45] "Engine< 1366"   "Engine>=1366"   "Seats=4,7"      "Power>=99.41"  
## [49] "Fuel_Type=Disl" "Fuel_Type=Ptrl" "Power< 99.41"   "Seats=5,8"     
## [53] "Engine< 1496"   "Engine>=1496"   "Engine< 1496"   "Engine>=1496"  
## [57] "Engine< 1354"   "Seats=6"        "Seats=4,5,7,8"

Los errores que se cometen en cada hoja y el error total del árbol:

printcp(arbolDecision)
## 
## Classification tree:
## rpart(formula = data[, 15] ~ ., data = data[, c(2:10, 13, 14)], 
##     method = "class", parms = list(split = "information"), cp = 0.004)
## 
## Variables actually used in tree construction:
## [1] Engine    Fuel_Type kmpl      Power     Seats    
## 
## Root node error: 1353/4711 = 0.2872
## 
## n= 4711 
## 
##           CP nsplit rel error  xerror      xstd
## 1  0.7139690      0  1.000000 1.00000 0.0229528
## 2  0.0177384      1  0.286031 0.28603 0.0139298
## 3  0.0162602      3  0.250554 0.27494 0.0136808
## 4  0.0110865      4  0.234294 0.24908 0.0130737
## 5  0.0096083      6  0.212121 0.21951 0.0123293
## 6  0.0088692     10  0.173688 0.19956 0.0117915
## 7  0.0081301     12  0.155950 0.19586 0.0116884
## 8  0.0066519     13  0.147820 0.18551 0.0113933
## 9  0.0059128     14  0.141168 0.17221 0.0109993
## 10 0.0051737     17  0.123429 0.16482 0.0107727
## 11 0.0048780     18  0.118256 0.15743 0.0105401
## 12 0.0044346     23  0.093865 0.14043 0.0099802
## 13 0.0040000     29  0.064302 0.12712 0.0095146
rpart.plot(arbolDecision, uniform=T);

Vamos a comprobar la precisión del árbol con los mismos datos con los que ha sido construido.

predArbolInfo1 = predict(arbolDecision, data[,c(2:10,13,14)])
# Redondeamos la predicción para obtener un resultado 0-1
predArbolInfo1 = round(predArbolInfo1)
tArb1 = table(predArbolInfo1[,1], data[,15])
tArb2 = table(predArbolInfo1[,2], data[,15])
medidaPrecision(tArb2)
## [1] 0.9745568

Obtenemos una precisión del 97% con los datos de entrenamiento. Ahora, vamos a ver como clasifica nuestro árbol a nuevos coches. Para ello usaremos el conjunto de datos de testing y calcularemos la medida de precisión del resultado.

predArbolInfo2 = predict(arbolDecision,dataTest[,c(2:10,14,15)])
# Redondeamos la predicción para obtener un resultado 0-1
predArbolInfo2 = round(predArbolInfo2)
tArb3 = table(predArbolInfo2[,1], dataTest[,13])
tArb4 = table(predArbolInfo2[,2], dataTest[,13])
#medidaPrecision(tArb3)
medidaPrecision(tArb4)
## [1] 0.9553525

Con el conjunto de datos de testing se obtiene una precisión del 95%, lo que nos indica que nuestro árbol de decisión realiza buenas predicciones y no sobreajusta el modelo.

Por otro lado, vamos a ver cuál es la precisión que obtenemos utilizando el criterio de Ginni. Construimos un nuevo árbol:

arbolDecisionGini <- rpart(data[,15] ~., data=data[,c(2:10,13,14)], cp=0.004, parms=list(split="gini"))
arbolDecisionGini
## n= 4711 
## 
## node), split, n, loss, yval, (yprob)
##       * denotes terminal node
## 
##   1) root 4711 1353 Baja-media (0.287200170 0.712799830)  
##     2) Engine>=1690 1548  291 Alta (0.812015504 0.187984496)  
##       4) Power>=164.2 860   65 Alta (0.924418605 0.075581395)  
##         8) Fuel_Type=Diesel 726   31 Alta (0.957300275 0.042699725) *
##         9) Fuel_Type=Petrol 134   34 Alta (0.746268657 0.253731343)  
##          18) Power>=180.515 96    9 Alta (0.906250000 0.093750000) *
##          19) Power< 180.515 38   13 Baja-media (0.342105263 0.657894737)  
##            38) Power< 177.235 7    0 Alta (1.000000000 0.000000000) *
##            39) Power>=177.235 31    6 Baja-media (0.193548387 0.806451613) *
##       5) Power< 164.2 688  226 Alta (0.671511628 0.328488372)  
##        10) Engine>=2117.5 420   68 Alta (0.838095238 0.161904762)  
##          20) Power< 151 379   37 Alta (0.902374670 0.097625330)  
##            40) Fuel_Type=Diesel 371   31 Alta (0.916442049 0.083557951)  
##              80) Engine< 2498.5 332   17 Alta (0.948795181 0.051204819)  
##               160) Power< 137 198    0 Alta (1.000000000 0.000000000) *
##               161) Power>=137 134   17 Alta (0.873134328 0.126865672)  
##                 322) Power>=138.555 118    1 Alta (0.991525424 0.008474576) *
##                 323) Power< 138.555 16    0 Baja-media (0.000000000 1.000000000) *
##              81) Engine>=2498.5 39   14 Alta (0.641025641 0.358974359)  
##               162) Engine>=2511 27    3 Alta (0.888888889 0.111111111) *
##               163) Engine< 2511 12    1 Baja-media (0.083333333 0.916666667) *
##            41) Fuel_Type=Petrol 8    2 Baja-media (0.250000000 0.750000000) *
##          21) Power>=151 41   10 Baja-media (0.243902439 0.756097561) *
##        11) Engine< 2117.5 268  110 Baja-media (0.410447761 0.589552239)  
##          22) Engine< 1796.5 22    0 Alta (1.000000000 0.000000000) *
##          23) Engine>=1796.5 246   88 Baja-media (0.357723577 0.642276423)  
##            46) Kilometers_Driven< 41183.5 43   11 Alta (0.744186047 0.255813953) *
##            47) Kilometers_Driven>=41183.5 203   56 Baja-media (0.275862069 0.724137931)  
##              94) kmpl< 16.755 137   52 Baja-media (0.379562044 0.620437956)  
##               188) kmpl>=14.09 72   26 Alta (0.638888889 0.361111111)  
##                 376) Power>=134.1 53   10 Alta (0.811320755 0.188679245)  
##                   752) Power< 147.555 45    2 Alta (0.955555556 0.044444444) *
##                   753) Power>=147.555 8    0 Baja-media (0.000000000 1.000000000) *
##                 377) Power< 134.1 19    3 Baja-media (0.157894737 0.842105263) *
##               189) kmpl< 14.09 65    6 Baja-media (0.092307692 0.907692308) *
##              95) kmpl>=16.755 66    4 Baja-media (0.060606061 0.939393939) *
##     3) Engine< 1690 3163   96 Baja-media (0.030350933 0.969649067)  
##       6) Seats=6 7    0 Alta (1.000000000 0.000000000) *
##       7) Seats=4,5,7,8 3156   89 Baja-media (0.028200253 0.971799747)  
##        14) Power>=132.16 7    2 Alta (0.714285714 0.285714286) *
##        15) Power< 132.16 3149   84 Baja-media (0.026675135 0.973324865)  
##          30) Engine>=1353.5 1243   75 Baja-media (0.060337892 0.939662108)  
##            60) Engine< 1366 29    0 Alta (1.000000000 0.000000000) *
##            61) Engine>=1366 1214   46 Baja-media (0.037891269 0.962108731) *
##          31) Engine< 1353.5 1906    9 Baja-media (0.004721931 0.995278069) *

Vamos a comprobar el error que comete el árbol de decisión y qué aspecto tiene:

labels(arbolDecisionGini, pretty=T)
##  [1] "root"                         "Engine>=1690"                
##  [3] "Power>=164.2"                 "Fuel_Type=Disl"              
##  [5] "Fuel_Type=Ptrl"               "Power>=180.5"                
##  [7] "Power< 180.5"                 "Power< 177.2"                
##  [9] "Power>=177.2"                 "Power< 164.2"                
## [11] "Engine>=2118"                 "Power< 151"                  
## [13] "Fuel_Type=Disl"               "Engine< 2498"                
## [15] "Power< 137"                   "Power>=137"                  
## [17] "Power>=138.6"                 "Power< 138.6"                
## [19] "Engine>=2498"                 "Engine>=2511"                
## [21] "Engine< 2511"                 "Fuel_Type=Ptrl"              
## [23] "Power>=151"                   "Engine< 2118"                
## [25] "Engine< 1796"                 "Engine>=1796"                
## [27] "Kilometers_Driven< 4.118e+04" "Kilometers_Driven>=4.118e+04"
## [29] "kmpl< 16.76"                  "kmpl>=14.09"                 
## [31] "Power>=134.1"                 "Power< 147.6"                
## [33] "Power>=147.6"                 "Power< 134.1"                
## [35] "kmpl< 14.09"                  "kmpl>=16.76"                 
## [37] "Engine< 1690"                 "Seats=6"                     
## [39] "Seats=4,5,7,8"                "Power>=132.2"                
## [41] "Power< 132.2"                 "Engine>=1354"                
## [43] "Engine< 1366"                 "Engine>=1366"                
## [45] "Engine< 1354"
printcp(arbolDecisionGini)
## 
## Classification tree:
## rpart(formula = data[, 15] ~ ., data = data[, c(2:10, 13, 14)], 
##     parms = list(split = "gini"), cp = 0.004)
## 
## Variables actually used in tree construction:
## [1] Engine            Fuel_Type         Kilometers_Driven kmpl             
## [5] Power             Seats            
## 
## Root node error: 1353/4711 = 0.2872
## 
## n= 4711 
## 
##          CP nsplit rel error  xerror      xstd
## 1 0.7139690      0   1.00000 1.00000 0.0229528
## 2 0.0177384      1   0.28603 0.29268 0.0140761
## 3 0.0162602      3   0.25055 0.27273 0.0136302
## 4 0.0155211      4   0.23429 0.27273 0.0136302
## 5 0.0073910      6   0.20325 0.22764 0.0125399
## 6 0.0072062      9   0.17886 0.17591 0.0111105
## 7 0.0059128     13   0.15004 0.16408 0.0107497
## 8 0.0044346     14   0.14412 0.15004 0.0103011
## 9 0.0040000     22   0.10791 0.13378 0.0097507
rpart.plot(arbolDecisionGini, uniform=T)

La precisión que se obtiene con los datos de entrenamiento es la siguiente:

predArbolGini1 = predict(arbolDecisionGini,data[,c(2:10,13,14)])
# Redondeamos la predicción para obtener un resultado 0-1
predArbolGini1 = round(predArbolGini1)
tArb5 = table(predArbolGini1[,1], data[,15])
tArb6 = table(predArbolGini1[,2], data[,15])
medidaPrecision(tArb6)
## [1] 0.9587796

Se obtiene un accuracy del 95%. Además, con el conjunto de datos de testing se obtiene la siguiente métrica:

predArbolGini2 = predict(arbolDecisionGini, dataTest[,c(2:10,14,15)])
predArbolGini2=round(predArbolGini2)
tArb7 = table(predArbolGini2[,1],dataTest[,13])
tArb8 = table(predArbolGini2[,2],dataTest[,13])
pArb1 = medidaPrecision(tArb7)
pArb2 = medidaPrecision(tArb8)
pArb2
## [1] 0.9329822

Comprobamos que tanto con el conjunto de entrenamiento como con el conjunto de prueba las precisiones que obtenemos son muy buenas también con Gini, en este caso se obtiene un 93% en test por lo que consideramos buenos árboles de decisión. Distinguiendo entre el árbol de información y el del criterio de Gini, seleccionamos el de información por cometer menos error y ajustarse más a las nuevas entradas que se le proporcionan.

Bosques aleatorios

Un bosque de árboles o bosque aleatorio es un conjunto de árboles de decisión tales que los nodos de cada árbol dependen de los valores de un subconjunto de variables muestreado aleatoriamente. Es decir, para construir cada árbol se escogen un subconjunto del total de las variables explicativas. Cuando tenemos el bosque construido lo que hacemos es “pasar” los nuevos datos por todos los árboles del bosque. Cada árbol clasificará al dato en una clase, y el bosque, finalmente, clasificará al dato en la clase en la que más árboles hayan coincidido.

Una vez que sabemos qué es un bosque aleatorio cabe preguntarse cuál es el número óptimo de árboles que debemos construir y cuántas variables queremos que se seleccionen para cada uno de los árboles. Para descubrir esto, probaremos con distintos valores de cada uno y nos quedaremos con el que mejor medida de precisión nos proporcione.

Durante todo el proceso de análisis se han evaluado cada una de las métricas. Al obtener una precisión del 97% en entrenamiento y 95% en testing en árboles de decisión, se ha hecho constatar que no seleccionaremos como mejor método para predecir el de bosques aleatorios o random forest ya que al obtener esa precisión tan alta, siempre será mejor tener un árbol que un bosque de árboles en cuanto a complejidad y coste. No obstante, se realizan algunos cálculos para estimar cual sería la precisión y las predicciones de este bosque.

set.seed(123)

modelo <- ranger(
  formula = as.factor(data[,15]) ~ .,
  data = data[,c(2:10,13,14)],
  num.trees = 10,
  seed = 123
)
predicciones <- predict(modelo, data = dataTest[,c(2:10,14,15)])
predicciones
## Ranger prediction
## 
## Type:                             Classification 
## Sample size:                      1173 
## Number of independent variables:  11
#test_rmse <- sqrt(mean((predicciones - as.factor(data[,13]))^2))
#paste("Error de test (rmse) del modelo (log): ", round(test_rmse,2))

Máquinas de vectores de soporte (SVM: Support Vector Machines)

Vamos a aplicar SVM a nuestros datos. Las SVM constituyen un método basado en aprendizaje para la resolución de problemas de clasificación y regresión. En ambos casos, esta resolución se basa en una primera fase de entrenamiento (donde se les informa con múltiples ejemplos ya resueltos, en forma de pares {problema, solución}) y una segunda fase de uso para la resolución de problemas. En ella, las SVM se convierten en una “caja negra” que proporciona una respuesta (salida) a un problema dado (entrada).

svmdata = data[,c(2:10,13,14,15)]

Lo primero que hacemos es separar el conjunto de datos en los conjuntos de train y test, estableciendo el 70% de ellos para train y el 30% para test.

ind <- sample(2,nrow(svmdata), replace= TRUE, prob = c(0.7,0.3))
trainSet <- svmdata[ind==2,]
testSet <- svmdata[ind==1,]

Realizamos el entrenamiento del modelo. Vamos a probar los distintos kernels para comprobar cual de ellos aprende mejor.

SVM con kernel radial

Suponiendo que una observación del conjunto de datos de test se encuentra alejada de una observación de entrenamiento en términos de distancia euclídea, el kernel radial tiene un comportamiento muy local, en el sentido de que sólo las observaciones de entrenamiento cercanas a una observación de test tendrán efecto sobre su clasificación. Es importante tener en cuenta que una mayor flexibilidad no tiene porque mejorar las predicciones debido a que un modelo muy flexible puede ajustarse demasiado a los datos de entrenamiento.

Primero, entrenamos el modelo:

model <- svm(Gama~., data = trainSet, kernel="radial")
prediccion <- predict(model, newdata= testSet[-12])

Ahora, se muestra el resultado mediante la matriz de confusión:

MC <- table(testSet[,12], prediccion)
MC
##             prediccion
##              Alta Baja-media
##   Alta        845         96
##   Baja-media  139       2222
accuracy <- (sum(diag(MC)))/(sum(MC))
accuracy
## [1] 0.928831

Como se puede observar, el modelo clasifica bien en un 92% de los casos, por lo que las predicciones son realmente buenas.

SVM con kernel lineal

El kernel lineal cuantifica la similitud de un par de observaciones usando la correlación de Pearson. Con un kernel lineal, el clasificador obtenido es equivalente a un support vector classifier. Se entrena y se muestran los resultados a continuación:

model2 <- svm(Gama~., data = trainSet, kernel="linear")
prediccion <- predict(model2, newdata= testSet[-12])

MC <- table(testSet[,12], prediccion) 
MC
##             prediccion
##              Alta Baja-media
##   Alta        848         93
##   Baja-media  140       2221
accuracy <- (sum(diag(MC)))/(sum(MC))
accuracy
## [1] 0.9294367

Se obtiene un accuracy del 92%, por lo que predice con bastante exactitud al igual que el kernel radial.

SVM con kernel sigmoidal

Ahora, con el kernel sigmoidal:

model3 <- svm(Gama~., data = trainSet, kernel="sigmoid")
prediccion <- predict(model3, newdata= testSet[-12])

MC <- table(testSet[,12], prediccion)
MC
##             prediccion
##              Alta Baja-media
##   Alta        848         93
##   Baja-media  142       2219
accuracy <- (sum(diag(MC)))/(sum(MC))
accuracy
## [1] 0.928831

Se obtiene un accuracy del 92%, por lo que predice con bastante exactitud al igual que el kernel radial y lineal.

SVM con kernel polinomial

El kernel polinómico de grado d (siendo d>1) permite un límite de decisión mucho más flexible. Cuando un support vector classifier se combina con un kernel no lineal, se obtiene un support vector machine. Se entrena el modelo y se muestran los resultados:

model4 <- svm(Gama~., data = trainSet, kernel="polynomial")
prediccion <- predict(model4, newdata= testSet[-12])

MC <- table(testSet[,12], prediccion)
MC
##             prediccion
##              Alta Baja-media
##   Alta        322        619
##   Baja-media   26       2335
accuracy <- (sum(diag(MC)))/(sum(MC))
accuracy
## [1] 0.8046638

Tanto el kernel radial como el lineal y el sigmoidal obtienen la misma métrica, por lo tanto esos hiperplanos se adaptan correctamente a nuestros datos.

Ajuste de los hiperparámetros del modelo

Muchos modelos, entre ellos los árboles de regresión, contienen parámetros que no pueden aprenderse a partir de los datos de entrenamiento y, por lo tanto, deben de ser establecidos por el analista. A estos se les conoce como hiperparámetros. Los resultados de un modelo pueden depender en gran medida del valor que tomen sus hiperparámetros, sin embargo, no se puede conocer de antemano cuál es el adecuado. Aunque con la práctica, los especialistas en machine learning ganan intuición sobre qué valores pueden funcionar mejor en cada problema, no hay reglas fijas. La forma más común de encontrar los valores óptimos es probando diferentes posibilidades.

  • Escoger un conjunto de valores para el o los hiperparámetros.

  • Para cada valor (combinación de valores si hay más de un hiperparámetro), entrenar el modelo y estimar su error mediante un método de validación

  • Finalmente, ajustar de nuevo el modelo, esta vez con todos los datos de entrenamiento y con los mejores hiperparámetros encontrados.

Regresión logística

El método glm de caret emplea la función glm() del paquete básico de R. Este algoritmo no tiene ningún hiperparámetro pero, para que efectúe una regresión logística, hay que indicar family = “binomial”.

registerDoMC(cores = 4)

Se configuran el número de repeticiones y las particiones para cada repetición:

particiones  <- 10
repeticiones <- 5

Ahora, se ajustan los hiperparámetros:

hiperparametros <- data.frame(parameter = "none")

set.seed(123)
seeds <- vector(mode = "list", length = (particiones * repeticiones) + 1)
for (i in 1:(particiones * repeticiones)) {
  seeds[[i]] <- sample.int(1000, nrow(hiperparametros))}
seeds[[(particiones * repeticiones) + 1]] <- sample.int(1000, 1)

Y por último, se define el entrenamiento y se ajusta el modelo:

# DEFINICIÓN DEL ENTRENAMIENTO
control_train <- trainControl(method = "repeatedcv", number = particiones,
                              repeats = repeticiones, seeds = seeds,
                              returnResamp = "final", verboseIter = FALSE,
                              allowParallel = TRUE)
# AJUSTE DEL MODELO
set.seed(34220)
modelo_logistic <- train(Gama ~ ., data = data[,c(2:10,13,14,15)],
                         method = "glm",
                         tuneGrid = hiperparametros,
                         metric = "Accuracy",
                         trControl = control_train,
                         family = "binomial")
modelo_logistic
## Generalized Linear Model 
## 
## 4711 samples
##   11 predictor
##    2 classes: 'Alta', 'Baja-media' 
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold, repeated 5 times) 
## Summary of sample sizes: 4240, 4239, 4240, 4240, 4239, 4241, ... 
## Resampling results:
## 
##   Accuracy   Kappa    
##   0.9266427  0.8214096

Empleando un modelo de regresión logística se consigue un accuracy promedio de validación del 93%.

k-NN

Ahora aplicamos la busqueda bayesiana para encontrar el mejor número de vecinos para el k-NN.

# Hiperparámetros
hiperparametros <- data.frame(k = c(1, 2, 5, 10, 15, 20, 30, 50))

set.seed(123)
seeds <- vector(mode = "list", length = (particiones * repeticiones) + 1)
for (i in 1:(particiones * repeticiones)) {
  seeds[[i]] <- sample.int(1000, nrow(hiperparametros)) 
}
seeds[[(particiones * repeticiones) + 1]] <- sample.int(1000, 1)
# DEFINICIÓN DEL ENTRENAMIENTO

control_train <- trainControl(method = "repeatedcv", number = particiones,
                              repeats = repeticiones, seeds = seeds,
                              returnResamp = "final", verboseIter = FALSE,
                              allowParallel = TRUE)
# AJUSTE DEL MODELO

set.seed(34220)
modelo_knn <- train(Gama~ ., data = data[,c(2:10,13,14,15)],
                    method = "knn",
                    tuneGrid = hiperparametros,
                    metric = "Accuracy",
                    trControl = control_train)
modelo_knn
## k-Nearest Neighbors 
## 
## 4711 samples
##   11 predictor
##    2 classes: 'Alta', 'Baja-media' 
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold, repeated 5 times) 
## Summary of sample sizes: 4240, 4239, 4240, 4240, 4239, 4241, ... 
## Resampling results across tuning parameters:
## 
##   k   Accuracy   Kappa    
##    1  0.8794749  0.7056852
##    2  0.8680960  0.6771780
##    5  0.8976852  0.7480572
##   10  0.8866888  0.7157664
##   15  0.8693664  0.6631850
##   20  0.8526423  0.6066279
##   30  0.8244122  0.5050079
##   50  0.7812349  0.3321058
## 
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was k = 5.
# REPRESENTACIÓN GRÁFICA
ggplot(modelo_knn, highlight = TRUE) +
  scale_x_continuous(breaks = hiperparametros$k) +
  labs(title = "Evolución del accuracy del modelo KNN", x = "K") +
  theme_bw()

Con un modelo k-NN con 5 vecinos se consigue un accuracy de validación promedio del 89’7%.

Árboles de decisión

Ahora, vamos a realizar el ajuste de hiperparámetros del árbol de decisión, que ya obtenía muy buenas predicciones.

# Hiperparámetros
hiperparametros <- data.frame(parameter = "none")

set.seed(123)
seeds <- vector(mode = "list", length = (particiones * repeticiones) + 1)
for (i in 1:(particiones * repeticiones)) {
  seeds[[i]] <- sample.int(1000, nrow(hiperparametros))
}
seeds[[(particiones * repeticiones) + 1]] <- sample.int(1000, 1)
# DEFINICIÓN DEL ENTRENAMIENTO
control_train <- trainControl(method = "repeatedcv", number = particiones,
                              repeats = repeticiones, seeds = seeds,
                              returnResamp = "final", verboseIter = FALSE,
                              allowParallel = TRUE)
# AJUSTE DEL MODELO
set.seed(34220)
modelo_C50Tree <- train(Gama~ ., data = data[,c(2:10,13,14,15)],
                    method = "C5.0Tree",
                    tuneGrid = hiperparametros,
                    metric = "Accuracy",
                    trControl = control_train)
modelo_C50Tree
## Single C5.0 Tree 
## 
## 4711 samples
##   11 predictor
##    2 classes: 'Alta', 'Baja-media' 
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold, repeated 5 times) 
## Summary of sample sizes: 4240, 4239, 4240, 4240, 4239, 4241, ... 
## Resampling results:
## 
##   Accuracy  Kappa    
##   0.981405  0.9544161
summary(modelo_C50Tree$finalModel)
## 
## Call:
## C50:::C5.0.default(x = x, y = y, weights = wts)
## 
## 
## C5.0 [Release 2.07 GPL Edition]      Mon Apr 18 11:41:54 2022
## -------------------------------
## 
## Class specified by attribute `outcome'
## 
## Read 4711 cases (31 attributes) from undefined.data
## 
## Decision tree:
## 
## Engine <= 1599:
## :...Seats5 <= 0:
## :   :...Power <= 99:
## :   :   :...Seats6 <= 0: Baja-media (128/1)
## :   :   :   Seats6 > 0: Alta (7)
## :   :   Power > 99:
## :   :   :...Fuel_TypeDiesel <= 0:
## :   :       :...Engine <= 1499: Baja-media (13)
## :   :       :   Engine > 1499: Alta (3)
## :   :       Fuel_TypeDiesel > 0:
## :   :       :...Seats8 <= 0: Alta (15)
## :   :           Seats8 > 0: Baja-media (2)
## :   Seats5 > 0:
## :   :...Engine <= 1343:
## :       :...kmpl > 18.33: Baja-media (1346/1)
## :       :   kmpl <= 18.33:
## :       :   :...kmpl <= 17.6: Baja-media (344)
## :       :       kmpl > 17.6:
## :       :       :...Power <= 77: Baja-media (60)
## :       :           Power > 77:
## :       :           :...Power <= 82: Alta (8)
## :       :               Power > 82: Baja-media (50)
## :       Engine > 1343:
## :       :...Engine <= 1364: Alta (30)
## :           Engine > 1364:
## :           :...Engine > 1496: Baja-media (745/4)
## :               Engine <= 1496:
## :               :...Engine > 1493:
## :                   :...Engine <= 1495: Baja-media (13)
## :                   :   Engine > 1495: Alta (18)
## :                   Engine <= 1493:
## :                   :...Fuel_TypeDiesel <= 0: Baja-media (77/5)
## :                       Fuel_TypeDiesel > 0:
## :                       :...Power > 66: Baja-media (276)
## :                           Power <= 66:
## :                           :...Power <= 64: Baja-media (24)
## :                               Power > 64: Alta (4)
## Engine > 1599:
## :...Power > 163.7:
##     :...Fuel_TypeDiesel <= 0:
##     :   :...Power > 180: Alta (96/9)
##     :   :   Power <= 180:
##     :   :   :...Power <= 177.01: Alta (7)
##     :   :       Power > 177.01:
##     :   :       :...Engine <= 1797: Alta (5)
##     :   :           Engine > 1797: Baja-media (26/1)
##     :   Fuel_TypeDiesel > 0:
##     :   :...kmpl <= 10.93:
##     :       :...kmpl <= 10.6: Alta (25)
##     :       :   kmpl > 10.6: Baja-media (10)
##     :       kmpl > 10.93:
##     :       :...Engine > 2200: Alta (301)
##     :           Engine <= 2200:
##     :           :...Engine > 2149:
##     :               :...Engine <= 2179: Alta (26)
##     :               :   Engine > 2179: Baja-media (14)
##     :               Engine <= 2149:
##     :               :...Power > 167.7: Alta (326/3)
##     :                   Power <= 167.7:
##     :                   :...Power <= 167.62: Alta (20)
##     :                       Power > 167.62: Baja-media (4)
##     Power <= 163.7:
##     :...Seats8 > 0: Alta (95)
##         Seats8 <= 0:
##         :...Power > 147.8:
##             :...kmpl > 18.7: Alta (7)
##             :   kmpl <= 18.7:
##             :   :...kmpl > 12.9:
##             :       :...Fuel_TypeDiesel <= 0: Baja-media (29)
##             :       :   Fuel_TypeDiesel > 0:
##             :       :   :...Engine <= 1984: Alta (3)
##             :       :       Engine > 1984: Baja-media (26/2)
##             :       kmpl <= 12.9:
##             :       :...Engine <= 2092:
##             :           :...TransmissionManual <= 0: Alta (15/3)
##             :           :   TransmissionManual > 0: Baja-media (2)
##             :           Engine > 2092:
##             :           :...Engine <= 2400: Baja-media (17)
##             :               Engine > 2400:
##             :               :...Engine <= 2773: Alta (6/1)
##             :                   Engine > 2773: Baja-media (6)
##             Power <= 147.8:
##             :...Engine <= 1997:
##                 :...Engine <= 1798: Alta (45/1)
##                 :   Engine > 1798:
##                 :   :...Power <= 139.01:
##                 :       :...Seats5 > 0: Baja-media (81/1)
##                 :       :   Seats5 <= 0:
##                 :       :   :...Year <= 2014: Baja-media (3)
##                 :       :       Year > 2014: Alta (3)
##                 :       Power > 139.01:
##                 :       :...kmpl <= 14.8: Baja-media (10)
##                 :           kmpl > 14.8:
##                 :           :...kmpl <= 16.8: Alta (32)
##                 :               kmpl > 16.8:
##                 :               :...Power <= 142: Baja-media (13)
##                 :                   Power > 142: Alta (6)
##                 Engine > 1997:
##                 :...Year <= 2011:
##                     :...Power <= 120.7: Alta (31/1)
##                     :   Power > 120.7:
##                     :   :...kmpl <= 14.49: Baja-media (23/1)
##                     :       kmpl > 14.49: Alta (4)
##                     Year > 2011:
##                     :...Power > 138.1: Alta (106)
##                         Power <= 138.1:
##                         :...Power > 136: Baja-media (9)
##                             Power <= 136:
##                             :...Power > 88.73: Alta (103/1)
##                                 Power <= 88.73:
##                                 :...Power <= 70.02: Alta (9)
##                                     Power > 70.02: Baja-media (4)
## 
## 
## Evaluation on training data (4711 cases):
## 
##      Decision Tree   
##    ----------------  
##    Size      Errors  
## 
##      57   35( 0.7%)   <<
## 
## 
##     (a)   (b)    <-classified as
##    ----  ----
##    1337    16    (a): class Alta
##      19  3339    (b): class Baja-media
## 
## 
##  Attribute usage:
## 
##  100.00% Engine
##   68.99% Seats5
##   58.01% kmpl
##   45.38% Power
##   28.27% Fuel_TypeDiesel
##   14.96% Seats8
##    6.26% Year
##    2.87% Seats6
##    0.36% TransmissionManual
## 
## 
## Time: 0.0 secs

Empleando como modelo un árbol simple C5.0, se consigue un accuracy promedio de validación del 98’1%.

Random Forest

El método ranger de caret emplea la función ranger() del paquete ranger. Este algoritmo tiene 3 hiperparámetros:

  • mtry: número predictores seleccionados aleatoriamente en cada árbol.

  • min.node.size: tamaño mínimo que tiene que tener un nodo para poder ser dividido.

  • splitrule: criterio de división.

Aunque caret también incluye el método rf con la función rf() del paquete randomForest, este último solo permite optimizar el hiperparámetro mtry.

# Hiperparámetros
hiperparametros <- expand.grid(mtry = c(3, 4, 5, 7),
                               min.node.size = c(2, 3, 4, 5, 10, 15, 20, 30),
                               splitrule = "gini")

set.seed(123)
seeds <- vector(mode = "list", length = (particiones * repeticiones) + 1)
for (i in 1:(particiones * repeticiones)) {
  seeds[[i]] <- sample.int(1000, nrow(hiperparametros))
}
seeds[[(particiones * repeticiones) + 1]] <- sample.int(1000, 1)
# DEFINICIÓN DEL ENTRENAMIENTO
control_train <- trainControl(method = "repeatedcv", number = particiones,
                              repeats = repeticiones, seeds = seeds,
                              returnResamp = "final", verboseIter = FALSE,
                              allowParallel = TRUE)
# AJUSTE DEL MODELO

set.seed(34220)
modelo_rf <- train(Gama ~ ., data = data[,c(2:10,13,14,15)],
                   method = "ranger",
                   tuneGrid = hiperparametros,
                   metric = "Accuracy",
                   trControl = control_train,
                   # Número de árboles ajustados
                   num.trees = 500)
modelo_rf
## Random Forest 
## 
## 4711 samples
##   11 predictor
##    2 classes: 'Alta', 'Baja-media' 
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold, repeated 5 times) 
## Summary of sample sizes: 4240, 4239, 4240, 4240, 4239, 4241, ... 
## Resampling results across tuning parameters:
## 
##   mtry  min.node.size  Accuracy   Kappa    
##   3      2             0.9507990  0.8806760
##   3      3             0.9509681  0.8810356
##   3      4             0.9509680  0.8810744
##   3      5             0.9509255  0.8810132
##   3     10             0.9499067  0.8785655
##   3     15             0.9486324  0.8754411
##   3     20             0.9471894  0.8719629
##   3     30             0.9449389  0.8663402
##   4      2             0.9646382  0.9138751
##   4      3             0.9637462  0.9117117
##   4      4             0.9634485  0.9109872
##   4      5             0.9633642  0.9108188
##   4     10             0.9611564  0.9054930
##   4     15             0.9588218  0.8998535
##   4     20             0.9565711  0.8945093
##   4     30             0.9516469  0.8825840
##   5      2             0.9750379  0.9390685
##   5      3             0.9743592  0.9374681
##   5      4             0.9741901  0.9370737
##   5      5             0.9736372  0.9357505
##   5     10             0.9719390  0.9316412
##   5     15             0.9698595  0.9266553
##   5     20             0.9668035  0.9192999
##   5     30             0.9612837  0.9060455
##   7      2             0.9835284  0.9597806
##   7      3             0.9829770  0.9584159
##   7      4             0.9830614  0.9586304
##   7      5             0.9826792  0.9577242
##   7     10             0.9816599  0.9552617
##   7     15             0.9796647  0.9504233
##   7     20             0.9780094  0.9464033
##   7     30             0.9733400  0.9349875
## 
## Tuning parameter 'splitrule' was held constant at a value of gini
## Accuracy was used to select the optimal model using the largest value.
## The final values used for the model were mtry = 7, splitrule = gini
##  and min.node.size = 2.
modelo_rf$finalModel
## Ranger result
## 
## Call:
##  ranger::ranger(dependent.variable.name = ".outcome", data = x,      mtry = min(param$mtry, ncol(x)), min.node.size = param$min.node.size,      splitrule = as.character(param$splitrule), write.forest = TRUE,      probability = classProbs, ...) 
## 
## Type:                             Classification 
## Number of trees:                  500 
## Sample size:                      4711 
## Number of independent variables:  30 
## Mtry:                             7 
## Target node size:                 2 
## Variable importance mode:         none 
## Splitrule:                        gini 
## OOB prediction error:             1.59 %
# REPRESENTACIÓN GRÁFICA
ggplot(modelo_rf, highlight = TRUE) +
  scale_x_continuous(breaks = 1:30) +
  labs(title = "Evolución del accuracy del modelo Random Forest") +
  guides(color = guide_legend(title = "mtry"),
         shape = guide_legend(title = "mtry")) +
  theme_bw()

Empleando un modelo random forest con mtry = 7, min.node.size = 2 y splitrule = “gini”, se consigue un accuracy promedio de validación del 98’4%.

SVM

El método svmRadial de caret emplea la función ksvm() del paquete kernlab. Este algoritmo tiene 2 hiperparámetros:

  • sigma: coeficiente del kernel radial.

  • C: penalización por violaciones del margen del hiperplano.

Se ajustan los hiperparámetros y se entrena el modelo con kernel radial:

# Hiperparámetros
hiperparametros <- expand.grid(sigma = c(0.001, 0.01, 0.1, 0.5, 1),
                               C = c(1 , 20, 50, 100, 200, 500, 700))

set.seed(123)
seeds <- vector(mode = "list", length = (particiones * repeticiones) + 1)
for (i in 1:(particiones * repeticiones)) {
  seeds[[i]] <- sample.int(1000, nrow(hiperparametros))
}
seeds[[(particiones * repeticiones) + 1]] <- sample.int(1000, 1)
# DEFINICIÓN DEL ENTRENAMIENTO
control_train <- trainControl(method = "repeatedcv", number = particiones,
                              repeats = repeticiones, seeds = seeds,
                              returnResamp = "final", verboseIter = FALSE,
                              allowParallel = TRUE)
# AJUSTE DEL MODELO

set.seed(342)
modelo_svmrad <- train(Gama ~ ., data = data[,c(2:10,13,14,15)],
                   method = "svmRadial",
                   tuneGrid = hiperparametros,
                   metric = "Accuracy",
                   trControl = control_train)
modelo_svmrad
## Support Vector Machines with Radial Basis Function Kernel 
## 
## 4711 samples
##   11 predictor
##    2 classes: 'Alta', 'Baja-media' 
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold, repeated 5 times) 
## Summary of sample sizes: 4240, 4240, 4239, 4240, 4240, 4240, ... 
## Resampling results across tuning parameters:
## 
##   sigma  C    Accuracy   Kappa    
##   0.001    1  0.9271074  0.8244331
##   0.001   20  0.9269789  0.8236833
##   0.001   50  0.9277433  0.8254514
##   0.001  100  0.9279981  0.8257346
##   0.001  200  0.9292287  0.8283934
##   0.001  500  0.9296964  0.8290602
##   0.001  700  0.9297806  0.8289960
##   0.010    1  0.9296963  0.8298518
##   0.010   20  0.9300766  0.8287067
##   0.010   50  0.9293980  0.8270251
##   0.010  100  0.9294410  0.8270782
##   0.010  200  0.9284219  0.8245447
##   0.010  500  0.9269788  0.8206876
##   0.010  700  0.9271481  0.8208105
##   0.100    1  0.9260447  0.8192948
##   0.100   20  0.9201869  0.8035070
##   0.100   50  0.9173425  0.7958958
##   0.100  100  0.9155593  0.7916444
##   0.100  200  0.9123322  0.7840549
##   0.100  500  0.9063028  0.7695067
##   0.100  700  0.9036705  0.7629402
##   0.500    1  0.9010397  0.7500010
##   0.500   20  0.8968771  0.7400260
##   0.500   50  0.8900000  0.7235634
##   0.500  100  0.8852446  0.7115387
##   0.500  200  0.8838005  0.7079148
##   0.500  500  0.8812541  0.7010558
##   0.500  700  0.8806599  0.6992638
##   1.000    1  0.8944579  0.7267584
##   1.000   20  0.8871119  0.7111117
##   1.000   50  0.8814231  0.6968666
##   1.000  100  0.8800224  0.6931900
##   1.000  200  0.8783665  0.6888309
##   1.000  500  0.8773057  0.6860360
##   1.000  700  0.8773906  0.6861693
## 
## Accuracy was used to select the optimal model using the largest value.
## The final values used for the model were sigma = 0.01 and C = 20.
modelo_svmrad$finalModel
## Support Vector Machine object of class "ksvm" 
## 
## SV type: C-svc  (classification) 
##  parameter : cost C = 20 
## 
## Gaussian Radial Basis kernel function. 
##  Hyperparameter : sigma =  0.01 
## 
## Number of Support Vectors : 878 
## 
## Objective Function Value : -13539.15 
## Training error : 0.057737
# REPRESENTACIÓN GRÁFICA
ggplot(modelo_svmrad, highlight = TRUE) +
  labs(title = "Evolución del accuracy del modelo SVM Radial") +
  theme_bw()

Empleando un modelo SVM Radial con sigma = 0.01 y C = 20, se consigue un accuracy promedio de validación del 93%.

Evaluación y comparación de modelos

Para poder determinar si un método es superior a otro, no es suficiente con comparar los mínimos (o máximos dependiendo de la métrica) que ha conseguido cada uno, sino que hay que tener en cuenta sus varianzas para determinar si existen evidencias suficientes de superioridad.

Al tratarse de modelos entrenados y validados sobre los mismos datos, mismas particiones y en el mismo orden (siempre que se haya asegurado la reproducibilidad mediante semillas), se pueden emplear métodos estadísticos para datos dependientes.

modelos <- list(KNN = modelo_knn, logistic = modelo_logistic,
                arbol = modelo_C50Tree, rf = modelo_rf,
                SVMradial = modelo_svmrad)

resultados_resamples <- resamples(modelos)

Se modifica el conjunto de datos para separar el nombre de cada modelo y las distintas métricas en la tabla:

metricas_resamples <- resultados_resamples$values %>%
                         gather(key = "modelo", value = "valor", -Resample) %>%
                         separate(col = "modelo", into = c("modelo", "metrica"),
                                  sep = "~", remove = TRUE)
metricas_resamples %>% 
  group_by(modelo, metrica) %>% 
  summarise(media = mean(valor)) %>%
  spread(key = metrica, value = media) %>%
  arrange(desc(Accuracy))
metricas_resamples %>%
  filter(metrica == "Accuracy") %>%
  group_by(modelo) %>% 
  summarise(media = mean(valor)) %>%
  ggplot(aes(x = reorder(modelo, media), y = media, label = round(media, 2))) +
    geom_segment(aes(x = reorder(modelo, media), y = 0,
                     xend = modelo, yend = media),
                     color = "grey50") +
    geom_point(size = 7, color = "firebrick") +
    geom_text(color = "white", size = 2.5) +
    scale_y_continuous(limits = c(0, 1)) +
    # Accuracy basal
    geom_hline(yintercept = 0.88, linetype = "dashed") +
    labs(title = "Validación: Accuracy medio repeated-CV",
         subtitle = "Modelos ordenados por media",
         x = "modelo") +
    coord_flip() +
    theme_bw()

metricas_resamples %>% filter(metrica == "Accuracy") %>%
  group_by(modelo) %>% 
  mutate(media = mean(valor)) %>%
  ungroup() %>%
  ggplot(aes(x = reorder(modelo, media), y = valor, color = modelo)) +
    geom_boxplot(alpha = 0.6, outlier.shape = NA) +
    geom_jitter(width = 0.1, alpha = 0.6) +
    scale_y_continuous(limits = c(0, 1)) +
    # Accuracy basal
    geom_hline(yintercept = 0.88, linetype = "dashed") +
    theme_bw() +
    labs(title = "Validación: Accuracy medio repeated-CV",
         subtitle = "Modelos ordenados por media") +
    coord_flip() +
    theme(legend.position = "none")

El modelo random forest consigue el accuracy promedio más alto, seguido muy de cerca por el árbol de decisión y SVM. Para determinar si las diferencias entre ellos son significativas, se recurre a test estadísticos. Como se comentó al principio, teniendo en cuenta las métricas obtenidas por cada uno de ellos, siempre es menos complejo y costoso seleccionar el árbol de decisión por delante de un bosque de árboles.

Test de Friedman para comparar el accuracy de los modelos

matriz_metricas <- metricas_resamples %>% filter(metrica == "Accuracy") %>%
                   spread(key = modelo, value = valor) %>%
                   dplyr::select(-Resample, -metrica) %>% as.matrix()
friedman.test(y = matriz_metricas)
## 
##  Friedman rank sum test
## 
## data:  matriz_metricas
## Friedman chi-squared = 177.79, df = 4, p-value < 2.2e-16

Para un nivel de significancia (α = 0.05), el test de Friedman sí encuentra evidencias para rechazar la hipótesis nula de que los clasificadores consiguen la misma precisión, sin embargo, no determina que par o pares son diferentes. Para identificarlos, se recurre a contrastes post HOC.

La función diff() del paquete caret recibe como argumento los resultados de validación de dos o más modelos extraídos con resample() y hace comparaciones por pares aplicando un t-test pareado con correcciones por comparaciones múltiples. Esta función no permite mucha flexibilidad en cuanto a las comparaciones, por lo que, una vez extraídos los datos con resample(), suele ser preferible emplear otras funciones disponibles en R.

# Comparaciones múltiples con un test suma de rangos de Wilcoxon
metricas_accuracy <- metricas_resamples %>% filter(metrica == "Accuracy")
comparaciones  <- pairwise.wilcox.test(x = metricas_accuracy$valor, 
                                        g = metricas_accuracy$modelo,
                                        paired = TRUE,
                                        p.adjust.method = "holm")

# Se almacenan los p_values en forma de dataframe
comparaciones <- comparaciones$p.value %>%
  as.data.frame() %>%
  rownames_to_column(var = "modeloA") %>%
  gather(key = "modeloB", value = "p_value", -modeloA) %>%
  na.omit() %>%
  arrange(modeloA) 

comparaciones

Acorde a las comparaciones por pares, no existen evidencias suficientes para considerar que la capacidad predictiva de los modelos es distinta.

Curva de ROC

Las curvas ROC (Receiver Operating Characteritic curve) permiten evaluar, en problemas de clasificación binaria, cómo varia la proporción de verdaderos positivos (sensibilidad o recall) y la de falsos positivos (1-especificidad) dependiendo del cutoff de probabilidad empleado en las asignaciones. El gráfico resultante es muy útil para identificar el cutoff que consigue un mejor equilibrio sensibilidad-especificidad. Además de esto, la curva ROC, en concreto el área bajo la curva (AUC), puede emplearse como métrica para evaluar modelos. Un modelo que clasifica perfectamente las dos clases tiene un 100% de sensibilidad y especificidad, por lo que el área bajo la curva es de 1. Un modelo que predice por debajo de lo esperado por azar, tiene un AUC menor de 0.5. Una condición necesaria para crear una curva ROC es disponer de la probabilidad de clases en las predicciones.

En caret, se puede sustituir la métrica accuracy empleada por defecto en problemas de clasificación y calcular en su lugar el AUC. Para ello, se tienen que indicar los argumentos summaryFunction = twoClassSummary y classProbs = TRUE en el control de entrenamiento. El segundo argumento es necesario porque el cálculo de la curva ROC requiere las probabilidades predichas para cada clase. Además del área bajo la curva, se calcula la sensibilidad y la especificidad para un cutoff de 0.5.

El paquete pROC contiene múltiples funciones para crear, representar y obtener métricas a partir de curvas ROC. Como argumentos se necesitan únicamente las probabilidades predichas para cada clase y la clase verdadera a la que pertenece cada observación.

# Se obtienen las probabilidades predichas para cada clase
predicciones <- predict(object = modelo_logistic,
                        newdata = dataTest[,c(2:10,13,14,15)],
                        type = "prob")
# Cálculo de la curva
curva_roc <- roc(response = dataTest$Gama, 
                 predictor = predicciones$Alta) 

# Gráfico de la curva
plot(curva_roc)

# Área bajo la curva AUC
auc(curva_roc)
## Area under the curve: 0.9471
# Intervalo de confianza de la curva
ci.auc(curva_roc, conf.level = 0.95)
## 95% CI: 0.9314-0.9629 (DeLong)

Con la regresión logística se consigue un área bajo la curva del 94%.

Ahora, veamos la curva ROC para el modelo k-nn.

# Se obtienen las probabilidades predichas para cada clase
predicciones <- predict(object = modelo_knn,
                        newdata = dataTest[,c(2:10,13,14,15)],
                        type = "prob")
# Cálculo de la curva
curva_roc <- roc(response = dataTest$Gama, 
                 predictor = predicciones$Alta) 

# Gráfico de la curva
plot(curva_roc)

# Área bajo la curva AUC
auc(curva_roc)
## Area under the curve: 0.9184
# Intervalo de confianza de la curva
ci.auc(curva_roc, conf.level = 0.95)
## 95% CI: 0.8995-0.9372 (DeLong)

Ahora, para el modelo de árbol de decisión:

# Se obtienen las probabilidades predichas para cada clase
predicciones <- predict(object = modelo_C50Tree,
                        newdata = dataTest[,c(2:10,13,14,15)],
                        type = "prob")
# Cálculo de la curva
curva_roc <- roc(response = dataTest$Gama, 
                 predictor = predicciones$Alta) 

# Gráfico de la curva
plot(curva_roc)

# Área bajo la curva AUC
auc(curva_roc)
## Area under the curve: 0.9939
# Intervalo de confianza de la curva
ci.auc(curva_roc, conf.level = 0.95)
## 95% CI: 0.9898-0.998 (DeLong)

El área bajo la curva del árbol de decisión es de 99%, por lo que se realiza un buen diagnóstico con este modelo.

GAM

En los modelos GAM se pueden aplicar funciones (lineales o no lineales) a cada uno de los predictores. El paquete gamlss implementa varias funciones smooth que suelen dar buenos resultados.

modelo_gam <- gam(Gama ~ s(Engine) + s(Power) + Transmission,
                family=binomial,
                data    = data)

summary(modelo_gam)
## 
## Call: gam(formula = Gama ~ s(Engine) + s(Power) + Transmission, family = binomial, 
##     data = data)
## Deviance Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.9569 -0.2039  0.1593  0.3037  2.5937 
## 
## (Dispersion Parameter for binomial family taken to be 1)
## 
##     Null Deviance: 5649.673 on 4710 degrees of freedom
## Residual Deviance: 2071.113 on 4701 degrees of freedom
## AIC: 2091.113 
## 
## Number of Local Scoring Iterations: NA 
## 
## Anova for Parametric Effects
##                Df Sum Sq Mean Sq  F value    Pr(>F)    
## s(Engine)       1 1089.5 1089.48 1254.821 < 2.2e-16 ***
## s(Power)        1   20.6   20.64   23.772  1.12e-06 ***
## Transmission    1    9.5    9.52   10.962  0.000937 ***
## Residuals    4701 4081.6    0.87                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Anova for Nonparametric Effects
##              Npar Df Npar Chisq    P(Chi)    
## (Intercept)                                  
## s(Engine)          3     351.69 < 2.2e-16 ***
## s(Power)           3     106.37 < 2.2e-16 ***
## Transmission                                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Debido a la incorporación de las funciones no lineales (smoothers), hay que ser cauto a la hora de interpretar los coeficientes y los errores mostrados en el summary. Los coeficientes y errores de los predictores con funciones smooth (Engine y Power) contemplan únicamente la parte lineal, ignorando la no lineal. En el caso de los predictores lineales (Transmission) sus errores se estiman asumiendo que los términos con smoother son fijos, y no contemplan la incertidumbre introducida al ajustar cada una de las funciones smooth.

Para conocer la contribución total (lineal y no lineal) de los predictores transformados por funciones smooth se puede emplear la función anova. Esta función calcula el impacto que tiene en el modelo (en términos de grados de libertad totales, AIC y significancia estadística) el eliminar cada uno de los predictores de forma secuencial.

modelo_gam
## Call:
## gam(formula = Gama ~ s(Engine) + s(Power) + Transmission, family = binomial, 
##     data = data)
## 
## Degrees of Freedom: 4710 total; 4701 Residual
## Residual Deviance: 2071.113
anova(modelo_gam, test= "F")
plot(modelo_gam)

AIC(logit3, modelo_gam)

El modelo GAM ha conseguido reducir todavía más el valor AIC y la distribución de sus residuos ha mejorado ligeramente.